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ABSTRACT 


The problem of potential flow around two-dimensional 
airfoils is solved by using a new singular integral method. The 
potential flow equations for incompressible potential flow are 
written in a singular integral equation. This equation is solved at 
N collocation points on the airfoil surface. A unique feature of 
this method is that the airfoil geometry is specified as an 
independant variable in the exact integral equation. 

Compared to other numerical methods, the present 
calculation procedure is much simpler and gives remarkable 
accuracy for many body shapes. An advantage of the present 
method is that it allows the inverse design calculation and the 
results are extremely accurate. Compared to other previous 
calculations, the present design solution is simpler, more accurate 
and does not use an iteration procedure. 
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CHAPTER 1 

REVIEW OF EXISTING METHODS 


1.1 Introduction 

A potential flow is one which is inviscid and 
irrotational. The irrotational condition implies that the 
velocity can be defined in terms of a potential function by: 

U=V<t> (1.1) 

When the problem involves a prescribed free stream flow over 
an arbitrary body, the velocity is commonly expressed as: 

WL+q, (1.2) 

where is the onset flow present when the body is not 

present and q is the disturbance velocity. In most cases U,*, is 

a uniform flow defined as parallel to the x-axis. When the flow 
is potential and incompressible, the Navier-Stokes equations 
reduce to the following equation: 

V 2 <|>=0 (1.3) 

The flow field is completely determined by kinematics, 
when the appropriate boundary conditions are specified. On the 
surface of the airfoil, the vector velocity is tangent to the 
surface and the disturbance velocity vanishes as the distance 
from the airfoil increases to infinity. 

When the flow is not symmetrical, we need an additional 
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condition which is given by the Kutta requirement. The Kutta 
condition states that the flow cannot go around the sharp 
trailing edge, but must leave the airfoil so that the upper and 
lower streams join smoothly at the trailing edge. This 
condition determines a unique value of the circulation. 

There are many techniques for calculating the 
incompressible potential flow around two-dimensional bodies; 
this chapter reviews several methods which have some 
similarity to the current singular integral method. 

1.2 Thin-Airfoil Theory 

The thin-airfoil theory uses several approximations in 
order to calculate the surface pressure distribution. The 
method of calculation is convenient for a rapid estimation of 
the velocity or pressure distribution over the airfoil. This 
theory had its beginnings in the early days of Thermodynamics 
with Munk [1], Birnbaum [2] and Glauert [3]. 

We assume that the airfoil is thin and that the camber 
and angle of attack are small. This suggests that the 
disturbance velocity is small compared to the free steam 
velocity . Since at the stagnation point this statement is 

evidently not true, the calculation is useful in regions which 
exclude this point. 

The differential equation governing the flow field is: 
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V 2 <>= 0 _ (1.4) 

The boundary condition along the mean camber line given by, 

»! = U M-u a (1.5) 

dy 00 d« 00 

where ti is the equation of the camber line, is the free 

stream velocity and a is the angle-of-attack. The far stream 
condition can be stated as: 

at infinity. (1.6) 

d» dy 

The solution for the flow field is obtained by superposing 
three problems: 

Problem 1 represents a thin symmetrical airfoil at zero angle 
of attack. 

Problem 2 represents the steady flow past a cambered airfoil 
of zero thickness at zero angle of attack 
Problem 3 represents a flat plate airfoil at an angle of attack. 

We first consider the problem of a thin symmetrical 
airfoil at zero angle of attack. The effect of thickness can be 
represented by a continuous distribution of sources along the 
x-axis. The disturbance potential can then be expressed by the 
following integral equation: 

i — 

Jqfl;) ,0 9 [(k-^) + g 2 ] 

0 


(1.7) 
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The proper source distribution denoted by q( % ) is determined 
by the surface boundary condition which gives: 

q ( h )=2 U ^L 
00 dK 

For the problem of the cambered airfoil of zero thickness 
we use a suitable distribution of vortices along the X-axis of 
the airfoil. The disturbance potential is given at the field point 
(x , y) by the integral relation: 

I r -1 



( 1 - 8 ) 


(1.9) 


Again the boundary conditions allow us to evaluate the axial 
distribution of vorticity. 

The flow field over a flat plate airfoil is solved by using 
a vortex distribution. The suitable distribution of vorticity is 

given as a solution of the following integral equation: 

1 

— f y(f) = -U a (1.10) 

2TT J s h-S 

After some mathematical manipulation, we obtain the vortex 
distribution: 

y ( 0 ) =2 U a 1-1:089 
00 sine 


The relationship between h and 0 is given by: 

H= — ( 1 + COS0 ) 

2 


( 1 - 12 ) 


The general solution for the flow over a thin airfoil is finally 
obtain by superposing the three previous calculations. This 
developments come from reference [4] 

1.3 Surface Distributions for Potential Flow 

The principle of this method is to sum sources, sinks, 
vortices or dipoles on the surface of the airfoil to form a flow 
field that satisfies the boundary conditions. The surface of the 
airfoil is approximated by N elements or panels with N points 
at which the singularity is to be evaluated. Each singularity is 
superposed with the uniform free stream and the resulting flow 
velocity must be tangent to each N elements of the airfoil at 
the points where the singularity is to be evaluated. All surface 
singularity methods, sometimes called panel methods, use the 
zero normal velocity on the surface to derive an integral 
equation for the singularity distribution. The evaluation of the 
proper distribution basically solves the problem and allows the 
computation of the pressure distribution. The most 
straightforward form to formulate this method is to use 
Green's theorem. The potential at any point P exterior to the 
airfoil can be expressed as: 
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<(>( P ) =-— ff 1 -^Mq) ds + -j-fJ <t>(q)— / J \ ds 

4ttJJ c r(p,q) dn « dn „ [r(p,q)J 


( 1 . 13 ) 


n denotes the normal to the surface at point q. The potential 
at a point p on the surface is given by: 


d 

1 

an 

q 

,r(p,q). 


ds 


( 1 . 14 ) 


Since d<))/dnq is prescribed, this is an integral for <j>(r). This 

equation represents a Fredholm integral equation, whose Kernel 
is given by: 


d 

1 

dn 

q 

r (p,q). 


( 1 . 15 ) 


However, the formulation that is more convenient is given by a 
surface distribution of unknown source strength 


(p) -n 


i 


(q) ds 


s r (P,q) 

This type of distribution gives a unique solution for the 
potential flow. Applying the boundary condition, we obtain the 
following integral equation: 


( 1 . 16 ) 


2tt g (p) - |J -^~T~ 

J { anplrlp, 


q) 


c (p) ds = - n . U 

P oo 


( 1 . 17 ) 



This is a Fredholm integral equation of the second kind whose 
kernel is: 
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(1.18) 


These equations are the basic formulation of the surface source 
density method to solve the problem of potential flow. The 
accuracy of this method is determined by the number of the 
elements used to approximate the body surface. For usual 
shapes, such as a two dimensional airfoil, 30 to 60 elements is 
sufficient. Usually the only interest is in the calculation of the 
surface velocity. The potential off the airfoil needs not to be 
calculated. 


Lift is deduced by means of a vorticity distribution on 
the surface. A conventional airfoil has a sharp trailing edge 
therefore for each angle of attack, there is a unique circulation 
that makes the potential flow velocity finite at the trailing 
edge. This condition is known as the Kutta condition. One 
technique is to put a vortex surface distribution on the surface 
of the airfoil. This method was proven to give the most 
accurate solution. The votex distribution would take the same 


N elements used for the source distribution and all the 


components will be summed up over the elements. The 
variation of the strength of the vortices is arbitrary, however a 
constant strength gives the most accurate solution. Another 
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technique is to place the vorticity distribution on the airfoil 
mean line. In this case only one vortex singularity needs to be 
placed on the mean line. This formulation has been reviewed by 
Maskew and Woodward [5] 

This method of solving the potential flow using a 
surface distribution is general and can be applied to any kind of 
bodies (two- dimensional, axisymmetric and even three 
dimensional shapes). Because of its versatility and accuracy, it 
has become the most popular technique for computing potential 
flows. This developments come from references [6] and [7]. 

1.4 Conformal Transformations 

Conformal transformations solve the potential flow field 
by using a complex transformed plane. The method simplifies 
the calculation of the flow field by solving for the flow field 
over a circle in the transformed plane which corresponds to the 
flow over a complicated airfoil shape in the real plane. 

Laplace's equation in the real plane transforms into Laplace's 
equation in the virtual plane and also the boundary conditions 
remain the same in both planes. The transformation maps 
points from the real plane using complex variables: 

Z = H + iV (1.19) 


into points on a transformed plane: 


(1.20) 
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The mapping is given by a function of the type: 

Z = £ + rj /£ + r 2 /£ 2 + r 3 /£ 3 + ... 

where r- { are real constants. The problem in the transformed 

plane reduces to the flow field over a cylinder. The exact 
solution is found by using a doublet, a uniform free stream and 
a value for the circulation which ensures the uniqueness of the 
solution. The solution is then transformed back into the real 
plane to give the pressure distribution on the airfoil. This 
method is limited to special airfoil profiles for which a 
conformal transformation exists. 

The most important conformal transformation is the 
Joukowsky transformation which leads to a family of airfoils 
known as Joukowsky airfoils. The Joukowsky transformation 
has the form: 

z - ; ♦ r/c 

and the inverse transformation used to transform back the 

solution into the real plane is given by: 

1 



In the transformed plane the center of the circle is displaced 
from the origin and the X-axis displacement is proportional to 
the thickness of the Joukowsky airfoil while the Y-axis 


( 1 . 21 ) 


( 1 . 22 ) 


( 1 . 23 ) 



displacement is determined by the camber. There are an 
infinite number of flows over a circle and the unique solution 
is found by invoking the Kutta condition which determined the 
position of the rear stagnation point on the circle. Another 
important aspect of the transformation is that at infinity the 
flows have exactly the same form and the angle of attack is 
also the same in either plane. More details about the Joukowsky 
transformation can be found in the reference [8]. Real airfoils 
are not Joukowsky airfoils, however the study of Joukowsky 
airfoil shapes give general trends for ideal flow over real 
airfoil shapes of similar thickness and camber. 



CHAPTER 2 

MATHEMATICAL FORMULATION 
AND NUMERICAL ALGORITHM 


2.1 Introduction 

The problem of predicting the pressure distribution about 
a two-dimensional airfoil has received considerable attention 
from various investigators. This study uses a recently 
developed [9] singular integral technique to solve the potential 
flow field over a two-dimensional airfoil. The calculation 
technique bears some resemblance to conventional singular 
integral methods as it reduces the formulation of the 
two-dimensional flow problem to the solution of an integral 
equation. However, this method has significant differences 
from other methods. It derives the integral equation by using a 
Fourier transform. Then, introducing a Taylor series expansion, 
the inverse transform is evaluated analytically. An important 
aspect is that the surface geometry of the airfoil appears 
explicitly in the integral equation. The advantage of this 
technique is that the formulation allows the inverse 
calculation to be easily performed, i.e. given a desired 
pressure distribution, the airfoil geometry can be found 
without iteration. 

Two main computer programs have been developed for 
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two-dimensional airfoils. The first program will compute a 
pressure distribution for arbitrary combinations of airfoil 
geometry and angle of attack while the second program will 
calculate an airfoil profile for a given pressure distribution. 

The programs were shown to be in good agreement with known 
results. The numerical results will be presented in a 
subsequent chapter. This chapter descibes the mathematical 
development and numerical solution associated with this new 
formulation of the potential flow problem. 

2.2 Formulation of the Equations 

Consider a symmetrical airfoil of moderate thickness at 
some angle of attack to a free stream. The outer flow is 
two-dimensional, inviscid and irrotational. As shown in Figure 
2.1 , the cartesian coordinate system is taken with the origin at 
the center point of the airfoil chord. The length of the airfoil 
chord is taken using non-dimensional variables to be equal to 2. 
The outer free stream at infinity is inclined at an angle of 

incidence a relative to the x-axis. The equation of the airfoil 
relative to the axis system is denoted by: 

y=n (x) (2.i) 

We introduce a disturbance velocity vector q («,y) due to the 
presence of the airfoil and we write: 
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U = U + q (x,y) (2.2) 

O© 

In the following analysis, represents a constant vector. In 

non-dimensional form, the magnitude of 11^ is equal to 1 . In a 

similar way, we introduce a disturbance potential for the 
disturbance velocity and we have: 

<t> = <l> 00 + <t> (2.3) 

It should be noted that the disturbance potential does not need 
to be small in the following formulation. 

The first step in the formulation is to divide the 
problem into an upper and lower plane as shown in Figure 2.2. 

The mathematical problem will then be solved independently 
for the lower and upper plane. In terms of the potential the 
flow field is described by the following mathematical 
formulation for the upper plane: 

Differential equation 

*JL + d 4> - o 
2 2 

dx dy 


(2.4) 
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Surface boundary condition 


dU I 


U(h) -oo < H < -1 
dH d« 

UJ(h) 1 < K < oo 


( 2 - 5 ) 


Far stream condition 

<)> u (x,y) equal (J)^ as the distance from the 
airfoil increases to infinity 

For the lower plane a similar set of equations applies and the 
potential function will be denoted by ^(H.y) . For symmetrical 
flows with a symmetrical airfoif only the upper half needs to 
be computed. However, in the general case, both the upper and 
lower plane are independently calculated. The parameters, 
U u (h) and ID u (h) are the y-components of the velocities along 
the centerline upstream and downstream of the airfoil 
respectively. For a symmetrical flow U u (h) and UJ u (h) are 
zero. For nonsymmetrical flows, we require that ( U u (k) , 
UJ u (h) ) and ( U*Ch) , UJ^h) ) match along the centerline, this 
ensures a unique solution and in effect replaces the usual Kutta 
condition. 

Finally we rewrite the mathematical formulation in 
terms of the disturbance potential and the disturbance upwash 
and downwash velocities defined by: 
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U(k) = sin a + v (k) 

0 <* 

LU(k) = sin a + w (») 


(2.6) 


The mathematical statement of the disturbance potential is 
given below for the upper plane: 

» 4> U ■ » 4> U . o 
2 2 
dH dy 


where, 


d± 


u 


L *y J 


= f U (H) 


y=rt(K) 


f U (») = 


(2.7) 


( 2 . 8 ) 


U(h) -oo < K < -1 

^ y 

( cosa + — ) -^3 sina 

dK d« 

UJ(h) 1 < H < oo 


-1 < H < 1 


u 


u 


and 0 as h 2 + y 2 -» oo 

dK dy 


(2.9) 


A similar set of equations exists for the lower plane. 


2.3 Solution bv Fourier Transform 

In order to find a solution for the above mathematical 
equations we use the Fourier Transform in the x-direction 
assuming that the function <|> satisfies the Dirichlet conditions 
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in every finite interval. The Fourier transform pair is given by: 

OO 

1 

<J> 

-OO 


1 f ~ -leu 

(s,y) = J <!>(«, y) e dH 

'V 2TT -OO 


OO 


^ if ISH 

4) (H,y) = _i— I o(s,y) e ds 

.oo 


The function <D(s,y) is the Fourier transform of <|>(H,y) . The 
differential equation for<D(s,y) is given by, 

-s 2 4» = 0 

2 

ay 


with the boundary conditions: 

[ do 

ay . 

y=ti(s) 


“ F(s) 


d<D 


oand _» o as y ->■ oo 


ay 


oo 


F(s) = J f(K) e 1 
v 2 it_oo 


-ISH 


dH 


f (K) = 


U(H) -oo < H < -1 

(cosa + -^ — ) -22 — sina — 1 < k < 1 

dH dH 
LU(h) 1 < h < oo 


(2.10) 

(2.11) 

( 2 . 12 ) 

( 2 . 13 ) 

( 2 . 14 ) 

( 2 . 15 ) 
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The solution is easily found to be: 
<t>(s,y) - 


Applying the inversion formula , the potential function is 
formally given by: 

-oo 


oo 


J 


1-00 


-IsKg-ri(s)) is(«-^) 


e 


d$ 




Although the solution is exact, the above equation is not 
useful from a computational perspective. However a simplified 
expression is found by expanding the exponential in a Taylor 
series over the transform variable s. A term by term 
integration then yields the following approximate expressions: 

OO 

|i.±Jr< 5 ) V «,MU )d 5 

-OO 

oo 

£ V K,H:5,d? 

3 -oo 


The Kernel functions of the above integral equations, accurate 
to O(e^) where e=t/c, are given by: 


K (H,y;$) 

n 


«l£ 

(H-$) 2 +(y-Ti ) 2 


K y (H,y;$) 


U~T) 

(H-£) 2 +(y-ri ) 2 


(2.16) 


(2.17) 


(2.18) 


(2.19) 


Although, (2.18) and (2.19) represent approximate solutions in 
general, it can be shown that they are exact on the surface. 

Also, it should be noted that if r|=0, and if the disturbance 
potential is set to zero in equation (2.15), we recover the usual 
thin airfoil equations. 

It may be observed that the equations (2.1 8) are singular 
Fredholm integral equations of the second kind and their 
Kernels are of difference type with a Cauchy singularity. For 

example, the behavior of K h is: K H (H,y; ± <*> , as 

The variation of f(£) in the Kernel function is continuous 
across the interval [-1 ,1 ] exept at the leading edge. The 
treatement of the singularity condition at requires the 
use of the Cauchy principal values. Another important aspect is 
that the Fredholm integral equation must be solved in an 
iterative manner since f(£) contains a term d<|>/dH which is 
unknown. The solution of the equations (2.18) by means of 
various optimized quadrature techniques is discussed in the 
following section. 

2.4 Outline of the Singular Integral Equation Solution 

E ms sslm s 

We wish to solve the integral equation (2.18) on the 
surface of the airfoil with y = ti(k). First we replace the 
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function f(H) by its respective values along the x-axis 
upstream and downstream of the airfoil and on the surface of 
the airfoil. The resulting integral equation then becomes: 



(2.20) 


The general procedure used in this calculation method is 
to represent the prescribed function f(«) along the airfoil 
surface by a linear combination of basis functions, in this case, 
powers of h: 

f(K)«2^a. K 1 -1 < h < 1 (2.21) 

NO 1 


or, 



rt 2 \ 1 

(1 -h ) a. k 

i 


-1 < K < 1 


( 2 . 22 ) 


Both functions have been tested, along with several other basis 
functions. The polynomial (2.21) yields better results and the 
remainder of this chapter will be devoted to computation 
methods using the polynomial (2.21). 

The coefficients aj are unknown and the function f(«) is 

a polynomial of degree n. We solve for the coefficients aj from 
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equation (2.20). These values will depend on the choice of the 
quadrature points Hj and the degree of the power series. The 

integral operand in the airfoil integral equation is not a 
continuous function on [-1,1] and the equation (2.20) is a 
singular integral equation which can be solved using the Cauchy 
principal values on the interval [-1,1]. Replacing f(«) by a 
power series, equation (2.20) may be written as: 



The last integral of the right hand side can be conveniently 
evaluated term by term by using the Cauchy principal values 
given in the appendix: 


1 



Finally any term of degree i can be written as a function of the 
term of degree i-1 : 


(2.23) 


(2.24) 
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In equation (2.23), the coefficients aj are unknown and 

we define a set of n+1 collocation points Kj as shown in figure 

2.3, where the integral equation (2.23) is to be evaluated at 
each point. The basic integral (2.23) can be expressed at any 
point Hj as a linear combination of the coefficients, aj . This 

coefficients are obtained by numerical integration and are a 
function of the geometry of the airfoil. Application of the 
above conditions gives a set of n+1 linear equations for the n+1 
unknown values of aj. 


2.4.1 Numerical Solution for Symmetrical Flow Field 

Consider first the case of a symmetrical airfoil at zero 
angle of attack. We see that the disturbance velocities U(h) 
and UI(h) are null. Inspection of (2.23) yields the expression: 


a? u _JL 


a» 


Tf 


1 

i=0 


a. I. 

i I 


where I j is the Cauchy integral of degree i, obtained from the 

integration of the basis functions. If we use the basis 
functions given by equation (2.21), and note that: 


(2.25) 


(2.26) 
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u 


f(H) = u +-£i-)^n 

dK dH 


we obtain: 



Solving of <|) K U , and substituting into equation (2.26) yields: 

h. i l 



j. 

if 


H 

JLDL 


= -1 


i=0 



(2.27) 


(2.28) 


(2.29) 


This equation applies at any Hj location between -1 and +1 . 
Applying this equation at n+1 points yields n+1 equations for 
the n+t unknowns.The coefficients of a j are called and 

define a matrix of dimension n+1 . The solution to the airfoil 
equations is finally obtained by inverting the matrix by 

Gaussian triangularization and the coefficients of the function 
f(w) are given by: 

aj = my -1 . (-l)j 

The computational time to solve the matrix by 
triangularization is less than 5 seconds using the CDC Dual 
Cyber System for a matrix of dimension no larger than 40. The 
computing time is independent of the geometry of the airfoil, 
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however it is dependent of the number of collocation points n. 
In order to reduce the round-off errors, we compute the matrix 
with double precision variables with 1 6 decimal accuracy 
which is useful for a matrix of dimension higher than 30 since 
the determinant of the matrix is a small number. The solution 
breaks down for n^40. This difficulty is caused by small 
valuesin the determinant when n is large. For high values of n, 
quantities of similar values are substracted in the calculation 
of the determinant which results in a dangerous loss of 
accuracy in the value of the determinant, consequently, the 
accuracy decreases as the number of quadrature points 
increases. On the other hand,' using a small number of points 
may not be sufficient to define the shape of the body especially 
very near the leading or stagnation point. However, for 
sufficiently small n the use of single precision variables 
allows a quicker computing time and the single precision 
calculation gives satisfactory results. Finally the pressure 
distribution is obtained by substituting the values a ; into the 


following expression: 


cp> 1 -[f(n)]‘ 


1 * 


1 
dx 


(2.30) 


with f(x) 


-1 

i=1 


a. x 

i 


i-1 


(2.31) 
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2.4.2 Numerical Procedure for Non-Svmmetrical Flow Field 
In addition to the solution of the basic potential flow 
problem over a symmetrical airfoil at zero angle of attack, the 
solution to the no n-sym metrical flow field has been 
incorporated into the numerical method. For this case, the 
upwash U(k) and downwash UJ(h) are unknown. Consequently, 
both the upper and lower planes must be solved simultaneously 
and the solutions must be matched along the cut. A possible 
methodology is to use some initial guess for U(») and UJ(h) and 
iterate until the change in the upwash and downwash is 
sufficiently small. The set of equations for this methodology 
are given below: 



An alternative approach, and one which proves to be 
superior, is to use an approximate representation for U(h) and 
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LU(h) and avoid the iteration procedure. This is accomplished 
by representing the actual airfoil by an "equivalent" Joukowsky 
airfoil of the same thickness for the purpose of obtaining U(h) 
and UJ(h) only. This procedure yields an approximate 
representation for the upwash and downwash. An extensive 
numerical investigation showed that the solution to (2.32) and 
(2.33) is sufficiently insensitive to this approximation to 
justify its application. A lengthy analysis is necessary to 
describe the flow about a Joukowsky airfoil and the details of 
the calculation are given in appendix A. The result for the 
upwash and downwash disturbance velocities are given by the 
following expressions: 



sina 


sina 


It should also be remembered that in the Joukowsky calculation 
the position of the leading edge is slighty different than -1 , 
for example it is equal to -1 .01 4405 for a thickness of 1 2%. 


(2.34) 


(2.35) 



The upwash and downwash integral terms of equations 
(2.32) and (2.33) are solved numerically by a simple 
trapezoidal rule. Boundary conditions on the airfoil surface are 
applied and the matrix elements are calculated using the same 
set of points distributed on the airfoil surface as for the 
symmetrical flow. Although there will be some additional 
terms in the calculation of the matrix elements, the numerical 
procedure remains the same. The matrix elements may be 
written as: 
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A similar set of matrix elements applies for the lower plane. 


2.5 Approximation of the Surface Boundary bv n Quadrature 
Points 

The integral equation described in the previous section is 
to be evaluated at a set of points Hj distributed on the airfoil 

surface as shown in figure 2.3. Special attention should be 
taken when choosing the quadrature points since the accuracy 
of the calculation is fully determined by the number and the 



distribution of the set of points. The quadrature of order n+1 
determines the number of unknown coefficients of the 
previously described function f(n). 

The spacing of the points must be small compared to the 
dimensions of the airfoil. In addition the local curvature of the 
airfoil should be considered in the point distribution. The 
proper distribution of the points over the airfoil surface will 
be largely a matter of experience and intuituion. As a first 
approach and one that proves to give satisfactory results, we 
use a set of equally spaced points along the x-axis of the 
airfoil. The first point is located at a distance d from the 
leading edge with succeding points spaced the same distance d. 
However two serious problems arise from the sharp corner at 
the trailing edge and from the large slope at the leading edge. 
These areas need to be defined by using a higher concentration 
of points. A higher order implementation which uses 
parabolically varying distances between points has been 
applied to the airfoil problem. A high concentration of points 
occurs at the leading and trailing edge and varies toward the 
central region of the airfoil where the distribution is sparce. 
However, for high order implementations, longer computing 
times will be required and loss of accuracy may occur from 
round-off errors in the matrix calculation. 



2.6 Solution of the Inverse Design Problem 
2.6.1 Introduction 

This section discusses an attempt to design by analytic 
means a class of airfoils using a similar methodology as for 
the direct problem. The design and development of aerodynamic 
bodies is usually an empirical procedure, based primarily upon 
the designer's experience and employing trial and error 
techniques. For the design problem, analytic solutions are not 
as developed as for the direct problem since it is more 
difficult and it involves the solution of a free boundary value 
problem. However, it is of great importance since for a 
desirable pressure distribution we can obtain the corresponding 
body shape. 

In the literature, solutions to the design problem are 
mostly based on iteration techniques due to the absence of 
exact mathematical solutions for free boundary value problems. 
Marshall [1 0] presented a technique that removes the free 
boundary element by a perturbation procedure. An analytical 
solution, using a surface source distribution, is obtained in the 
form of integral equations. Nevertheless, the calculation 
method uses an initial guess and involves an iteration 
procedure. Zedan and Dalton [1 1] presented a method which 
employs an axial source-sink distribution, with constant 
element strength, to obtain a solution to the design problem. 



The method proves to be accurate and converges, but it uses an 
iteration procedure and the method is also limited to bodies 
that do not present a sudden change in the slope of the meridian 
line. This present study does not requires any iteration and 
even less computational time is necessary than for the direct 
calculation. 

2.6.2 Mathematical Formulation 

In this section, the basic equations for the design 
problem with uniform flow field are derived. In this study, the 
method uses the surface velocity instead of the pressure as the 
prescribed distribution. The design problem can be stated as: 
given a surface velocity, what is the body shape that would 
produce this velocity distribution? 

Again consider an inviscid incompressible flow over a 
two-dimensional airfoil. Equation (2.20) of the previous section 
remains applicable since the flow conditions remain the same: 
This equation is repeated for convenience. 



We are now faced with the problem of finding the shape 
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prescribed by diyd£ given the surface velocity. As before, the 
function f(») may be expressed by a linear combination of 
powers of h. 

f(H)= A a. K 1 -1 < H < 1 

i=o 1 

where f(») is given by: 

f (h) = ( cosa + *24. )^2L - S jna = M. 

d« dK dy 

Equation (2.38) is evaluated at a set of n+1 quadrature points, 
which give a set of n+1 linear equations solved by Gaussian 
elimination, to obtain the n+1 values aj. It should be noted that 

the x-disturbance velocity is actually the unknown at this point 
since the airfoil shape is still unknown. This component is 
determined by inserting the coefficients aj in the integral 

equation (2.39). A similar procedure as the one used in the 
direct problem allows us to calculate the integrals of equation 
(2.37) by introducing the Cauchy principal values. 

In order to gain better accuracy, the x-disturbance 
velocity is evaluated at 200 points along the chord length and 
the slope on the airfoil surface is ultimately given by the 
equation: 


(2.38) 


(2.39) 
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The treatment of the inverse design problem has been 
restricted to uniform flow at a zero degree of angle of attack 
which requires a less sophisticated approach since the location 
of the stagnation point is known. 

Results of both the design and analysis problem 
determined by the procedure outlined above are presented in the 
next section. 



CHAPTER 3 

RESULTS FOR THE ANALYSIS MODE 


3.1 Introduction 

In this chapter, a series of numerical calculations for 
different airfoil geometries are presented. The results are 
validated by comparison to numerical solutions and analytical 
solutions when they exist. A code which was recently 
developed at NASA Langley [12] has been selected for purposes 
of verification of the present method. This code uses a 
spectral multigrid technique and has been extensively validated 
with finite difference schemes. In addition to this numerical 
verification, the present results are compared to analytical 
solutions for elliptical and Joukowski airfoils. 

Data are presented in terms of the pressure coefficient 
Cp, which is the quantity of usual aerodynamic interest. It is 
defined, in general, as: 

P-P^ 

Cp = — (3.1) 


where p denotes the local pressure. In incompressible 
potential flow it is related to the velocity U by : 


Cp = 1 



(3.2) 
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The formulation of the problem, presented in Chapter 2, has 
been tested for the flow over a 12% thick elliptic airfoil, a 
NACA 001 2 airfoil and a 1 2% thick Joukowski airfoil. A 
description of the results follows. 

3.2 Implementation of the Quadrature Points 

Flows have been computed using both equally spaced 
points and a higher order implementation. In the first method, 
the distribution of the points is simply determined by using a 
constant value A for the distance between two consecutive 
points. It should be also noted that the distance between the 
leading edge and the first collocation point as well as the 
distance between the last point and the trailing edge is equal 
to A. 

The second method uses a geometrically increasing grid. 
The following equation is applied to determine the spacing 
between two consecutive points: 

8k j = A rJ (3.3) 

The resulting distribution is shown in Figure 3.1 . Moving away 
from the leading edge, each one-dimensional grid spacing is 
made r times larger until the center of the airfoil is reached. 

The parameters A and r are constant values and j denotes the 
j th interval. In order to evaluate the two unknown A and r, we 
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CP = 1- F 


(1 + e)' 


1 - e 2 


H 


1-K 


(3.6) 


where e denotes the thickness ratio. 

Figure 3.2 and figure 3.3 compare the analytical solution 
with the calculated solution using 20 and 24 points distributed 
on the airfoil surface. Figure 3.2 uses the geometrically 
increasing grid and Figure 3.3 uses equally spaced points. The 
free stream velocity is parallel to the x-axis of the ellipse and 
as can be seen, the two plots are graphically indistinguishable 
for N=24. Positions of the points are shown in figures 3.4 and 
3.5. Figure 3.2 and 3.3 are representative of several other 
calculations which were made using using a larger and smaller 
number of points and various types of grid points distributions. 


3.4 NACA 0012 Airfoil 

The airfoil profile is given by the equation: 

y = 1.2( .29697k - .12600k - .35160k 2 +.2840h 3 
-.10150k 4 ) 1>k>0 (3.7) 

Two calculated pressure distributions are shown for the NACA 
0012 airfoil. One was calculated by using the methodology 
discussed in this study . The other distribution was obtained 
from a NASA computer code [8] and used as a comparison. As 
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POSITION OF THE QUADRATURE POINTS 
N- 24. 






illustrated in figures 3.6, 3.7 and 3.8, the two results are 
indistinguishable over the central region. Agreement with the 
NASA computer code can be improved at the trailing edge by 
adopting the variable grid scheme given by equation (3.3). This 
improvement can be seen by comparing figures 3.6 and 3.7 , but 
it should be noted that a slight loss of accuracy occurs in the 
region of the leading edge when using the geometrically 
increasing grid. Positions of the points are shown in figures 
3.9 and 3.10 for the distributions of figures 3.7 and 3.8. 

The calculations were repeated for different sets of 
quadrature points and the method has proven to give consistent 
results over a range of 1 5 to 30 points distributed over the 
airfoil surface. For N greater than approximatly 35, the 
accuracy begins to decrease due to the increasing matrix 
round-off error. A number of points smaller than 1 5 does not 
give an accurate description of the airfoil geometry. The 
calculated pressure coefficient exhibits a small repeated error 
very near the leading edge. This behavior can be explained by 
the difficulty that occurs when fitting a polynomial function 
over the region of large velocity gradients, e.g. the leading edge 
region. Slight changes in the locations and the number of the 
quadrature points can improve the accuracy of the curve. 
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Figure 3.9 
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Figure 3.10 



3.5 NACA 001 2 at 4 and 1 0 Degrees Angle of Attack 

Figures 3.1 1 , 3.12 and 3.13 show pressure distributions 
on the airfoil at 4 degrees angle of attack and figure 3.14 
shows calculated results for 10 degrees. Two curves are 
shown for each figure. One corresponds to the pressure 
distribution on the lower plane and the other is the pressure 
distribution on the upper plane. It should be remembered that 
the distributions for both the lower plane and the upper plane 
are independantly calculated and the solutions are matched 
along the x-axis upstream and downstream of the airfoil. 

In figure 3.1 1 , a equally spaced grid has been used while 
in figures 3.12 to 3.19 a geometrically increasing grid has been 
used. In figure 3.1 1 , the calculated distribution and the 
distribution obtained from the NASA code [9] are virtually 
identical. Again agreement with the NASA code is excellent in 
figures 3.1 2 and 3.13. Note that for a slight change of angle of 
attack from 0 to 4 degrees, the maximum peak of the pressure 
distribution experiences a change from approximately -0.5 to 
-1 .5. For 1 0 degrees angle of attack shown in figures 3. 1 4 and 
3.15, the upper plane calculation gives reasonably accurate 
results, while the lower plane calculation gives almost 
identical results compared to the data [13]. Positions of the 
collocation points are given by the variable grid scheme 
described in section 3.2 . The calculated Cp is slightly less 
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than the reference data in most of the upper central plane 
region. It should be noted that the trailing and leading edge 
regions are in close agreement with the data [9]. 

3.6 Joukowski Airfoil 

Figure 3.16 shows the computed pressure distribution for 
the case of the symmetrical Joukowski airfoil in a steady flow 
at a zero degree angle of attack. A comparison is made with 
the calculated pressure distribution obtained by using a 
transformed plane and the Joukowski transformation. Details 
about the calculation procedure can be found in reference [14]. 

The polynomial surface pressure distribution deviates 
quite seriously over the region of the negative high pressure 
peak, toward the leading edge. The correlation is quite 
reasonable on the surface of the right half airfoil plane toward 
the trailing edge and the general trend of the C p polynomial 

curve is in good agreement with the calculated curve. The main 
problem occurs in the region of the negative high pressure peak 
where the pressure distribution is overpredicted. Extensive 
calculations were made in an effort to improve the agreement. 
Different polynomial approximation schemes were used and the 
analytical solutions were carefully checked. In all cases, the 
present method consistently overpredicted the negative 
pressure peak. It has been concluded that this error is probably 



due to the error introduced in the conformal mapping which 
produces a slightly displaced leading edge. As noted in section 
2.4.2, this error is approximatly 1 .5% for the 12% thick 
Joukowski airfoil. It should be noted, that similar 
disagreements have been observed by other investigators [12]. 

Figures 3.17, 3.18 and 3.19 show pressure distribution 
calculations made for different grid point distributions. 

Positions of the points are shown in figures 3.20 and 3.21 . An 
important aspect is that the computation procedure gives 
consistent results using different grid point distributions. 

The essential difference between a cusped trailing edge 
and a trailing edge of finite angle is evident from a comparison 
of figure 3.7 (NACA 0012 airfoil) and figure 3.18 (Joukowski 
airfoil). The behavior of the flow at the trailing edge is 
accuratly calculated by the polynomial method as shown in 
figure 3.18 and the correlation for both the polynomial and 
theoretical curves agree well in the trailing edge region. 

The calculations for the flow about the symmetrical 
Joukowski airfoil of thickness 12% were repeatead at 4 and 1 0 
degrees of angle of attack. The calculated pressure 
distributions are compared with the analytic solution in 
figures 3.22 and 3.23. For the lower curve of figure 3.22, both 
pressure distributions are virtually identical. The largest 
disagreement occurs near the negative high pressure peak. 
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Nevertheless, it can be seen that the calculated pressure 
distribution is in close agreement with the theoritical 
distribution near the trailing edge. Slight changes in the 
locations of the collocation points do not significantly affect 
the shape of the polynomial curve. 

3.7 Inverse Design Results 

The problem of solving for the body shape given the 
surface velocity distribution uses essentially the same 
approach as the direct problem. The calculation procedure is 
described in detail in chapter 2. However, the technique is not 
based on an iteration technique as most design methods found 
in the literature. 

The surface velocity distributions used were exact 
solutions for the cylinder, the elliptic airfoil and the 
Joukowski airfoil and a numerical approximation for the NACA 
0012 airfoil. The calculated body shape is compared to the 
exact body to evaluate the accuracy of the method. 

Figures 3.24, 3.25, 3.26 and 3.27 show the Y-component 
of the surface velocity for the 4 described airfoils. Figures 
3.28, 3.29, 3.30 and 3.31 show the calculated and exact shape 
for the 4 airfoils. Using 24 equally spaced quadrature points, 
both the calculated and the exact shape are indistinguishable 


for the cylinder, the elliptic airfoil and the Joukowski airfoil. 

In the NACA 0012 airfoil case, the agreements for both curves 
are quite good. The calculated shape is slightly underestimated 
in the region of larger thickness. However, we should 
remember that the surface velocity distribution for the NACA 
0012 airfoil is not exact but obtained from a NASA code [12]. 
The imprecision in the NASA data could cause the small but not 
negligible errror of the calculated shape. Another 
consideration is that the error occurs in the high pressure peak 
region which also presented an error for the direct calculation. 

The design calculation presents less error everywhere as 
compared to the direct calculation and particularly near the 
leading and trailing edge. 
























CHAPTER 4 

CORRECTION FACTOR FOR THE 
COMPRESSIBLE CALCULATION 


4.1 Introduction 

This section presents an analysis for the problem of 
predicting the surface pressure distribution over a 
two-dimensional airfoil in a steady compressible potential 
flow. The flow field under consideration is inviscid, 
irrotational and the outer free stream velocity is limited to 
Mach numbers less than one. The main purpose of this section 
is to develop a numerical procedure that can be applied when 
the incompressible assumption is not valid. 

The solution procedure is similar to the incompressible 
flow calculation described in the preceding sections. The 
solution is reached in two steps. The initial pressure 
distribution is obtained by using the incompressible flow 
calculation and the solution is converted into the corresponding 
compressible solution by means of subsequent iterations which 
take into account the compressibility effect. The procedure is 
repeated until the solution converges. Details about the 
calculation method are fully described in the next sections. 
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4.2 Potential calculation method for compressible flow 
The equations are substantially modified to take into 
account the compressibility factor. However, the calculation 
procedure remains the same. The problem is divided into an 
upper and lower plane as shown in Figure 2.2 and the 
mathematical problem will be solved independantly for the 
lower and upper plane. The outer flow field is described by the 
following nonlinear system of equations for the upper plane: 

<|> u * H U (<|>) 

^ 1 - f(K) 

dy 

y=s(x) 


where, 
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<t> + 2 <t» <t> <t> 

hh x y xy 


+ $?.$ 


y yy 
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( 4 . 2 ) 


+ r±[ (<t> 2 + 4,2 _ ! ) ( <(, +(}) )] 
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On the solid surface of the airfoil the total velocity has to 

satisfy the tangency condition: 

f U(x) -go < « <-a 
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»y J 


y=s(x) 
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dd> as 
dx ax 


LU(x) 


-a < x < a 


a < x < +<» 


( 4 . 3 ) 


( 4 . 4 ) 


A similar set of equations can be derived for the lower half 
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plane. 

A successive approximation approach similar to the 
Rayleigh-Janzen method (2) will be adopted and we will iterate 
upon the solution for the compressible flow field. In order to 
apply this procedure, the maximum local Mach number will be 
restricted to values less than one. Using this approach and 
denoting each iteration with the superscript (n), equation 4.1 
and 4.2 can be rewritten as shown below: 

O < n+1 > 2 u (n) 

VV = M H (K,y) 


(4.5) 


d<t> 

»9 J 


y=$(«) 


U(h) -OO < K <-a 
-a < h < a 

dH dK 

W(H) a < k < +°° 


(4.6) 


The advantage of this approach is that for each 
approximation, the equations are linear and we can exploit 
several analytical techniques. For each approximation, we can 
decompose the potential function into the known freestream 
value plus the disturbance due to the airfoil. This follows from 
the linearity of the differential equation, and does not imply a 
small disturbance approximation. 
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We have: 


<t>(H,y)*<|> +<Mh, y) 

OQ 


U(k)* U + U (k) 
oo 


LU(k) = U) + LU (K) 
oo 


(4.7) 

(4.8) 
(4.9) 


Where U (h) and UJ (h) represent the unknown upstream 
and downstream influence of the airfoil. Substituting these 
expressions into equation 4.5 and 4.6 yields, 


V 2 <t> 


U <n * u , (n) 

u 2 u 

=M H (K,y) 


(4.10) 


u 


(n) 


a<j) 


«y J 


(n) 

= f U (h) 


(4.11) 


y=s(n) 


<j>H , <|>y -> 0 as | H I oo 


(4.12) 
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f U (H)= 


f u (x) is now 

given 

by: 

U(x) 


x<-a 

u 

J u 


( <|>x + $ ) 

ds 

- U 1 x l<a 

°°x 

dx 

oo 

W(X) 


x>a 


It will be advantageous to rewrite H(K,y) as shown below, 
H(x,y) = H m (K,y) + Hj(x,y) 

where, 

1 / , d A d , ..2 


H (x,y) = — ( $ — + <(> — )U 
m 2 x dx y dy 

(u 2 -n v 2 *] 

l)^= + <J>^ 

k y 


The system of equations given by 4.10 through 4.16 can 
now be solved in both the upper and lower planes. At each 
iteration, the solution is matched along the x-axis for |x|>a. 

This constraint, along with the requirement given below, 
uniquely determines the flow field, 

U(x) 0 as I x I -> <*> 


( 4 . 13 ) 


( 4 . 14 ) 

( 4 . 15 ) 

( 4 . 16 ) 


( 4 . 17 ) 


UJ (x) -¥ 0 as |x| oo 


( 4 . 18 ) 
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U) (+a) = 


dd) 

. ay . 


«=+a 


(4.19) 


4.3 Iteration procedure for the symmetrical flow case 

In order to illustrate the salient feature of the method, 
the less complicated case of a symmetrical flow will be 
considered. For this case, the upper and lower problems 
uncouple, i.e. , 


u (H) = ID (h) * 0 (4.20) 

and we can drop the superscript (u). The basic, n=0, 
approximation corresponding to the incompressible case is 
given in chapter 2. The incompressible flow solution serves as 
the basic, n=0, solution to the compressible flow problem. To 
compute a second approximation, the following system must be 
solved. 


2~* 2 * 2 (1) 

V 2 <t> = M H^'(H f y) 



,( 2 ), , 

= f (h) 


L ay J 


y=s(H) 


(4.21) 

(4.22) 


J2) J2) 

in i 5>y — >0 a$|Hl->oo (4.23) 


In equation 4.21 we used the result V 2 <t>^ ) = 0, to eliminate the 
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Hy^ ^(x,y) term. Taking the Fourier transform and eliminating 
the 0(e 2 ) term, we find the following solution. 

_(2) „(2) _ oo+oo _ 

d± 3 i!H__ + ML_r f H illi + iLl d ^ d 

dx dH 2H J J Mir r _ 1 


0 -oo 


(4.24) 

The real advantage of this approach is that the surface 
integral in equation 4.24 can be simplified to a line integral. 

This is accomplished as follows: Using the operator form of 


H m (£/n) , applying integration by parts, and invoking the 

Integral Mean Value Theorem to remove U 2 , the surface 
integral reduces to, 


oo+oo +oo 

D 2 J J V. K d ^ dr| + J ( 

0 -OO — oo 


a „ 2 ) 

dx dri y=s(K) 




In this expression, Kplx^^ri) is the kernel function given in 

equation 4.24 and K k = dKg/dx. The advantage of this 

re-formulation now becomes apparent. Applying Green's 
theorem to the surface integral, we see that it becomes zero, 
and the final result becomes, 
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Equation 4.25 represents a Fredholm integral equation which 
can be solved for d<)>/dH. The term in brackets represent the 
first compressibility correction to the basic (incompressible) 
flow solution. 

Successive approximations can now be determined 
immediately, once the differential system is stated. For 
example, the third approximation is given by, 


2 J3) 

V^<j> = M 


H (H,y) + H («,y) 
m T 


which becomes, 
J3) 


v< =m 2 h' 2 ' *M 4 £J-<u 2 -n (2) H (,) 

_ 2 - 


m 
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The associated surface boundary condition is 
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= f (h) 


y=s(«) 
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(4.27) 


(4.28) 
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and the solution is given by: 


J3) JD 
d<t> _ d<l) + 

dH dX 




K (H t y£) 

ri 




(rD u 




2 

-1) K (x,y;£) d£ 

K 


(4.29) 


In equation 4.29, K K (x,y,£) is determined by: 


K H (x,y£) = 


h-L 


(x-^) 2 +(y-s) 2 


(4.30) 


Higher approximation, i.e., n>4, can now be found by inspection, 
once the differencial equation is written. 


4.4 Numerical results 

A numerical result is presented for a symmetrical flow 
over the airfoil NACA 0012. The calculation is validated by 
comparison with another numerical computation from the 
computer code developed at NASA Langley [1 2]. 

An abbreviated iteration scheme was adapted for solving 
the compressible flow. This scheme is different from the 
derived relation given by equation (4.29) which was used in 
reference [9]. It includes only the first correction term given 
by equation (4.25), however, an iteration in the computer 
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program is performed on the integral equation until the 
solution reaches a converged value. The series of equations 
used for this approach are shown below for n<3. 


J2) J1) 

11 = 11 + J_f 

IT J 


dH dH 


^(3) Jl) 

11 =11 + _L 
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dH dH TT 
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(4.31) 


(4.32) 


(4.33) 


A subroutine was added to the code to compute the 
compressible flow from the results of the incompressible 
calculation. Results are presented in terms of the 
compressible pressure coefficient which is given by: 


C 

P 



1 + -*4-(U -q ) 
„ 2 00 
2a 

o© 


Y-l 


-1 


(4.34) 


Figure 4.1 shows the numerical solution using a 
polynomial of degree 20 with equally spaced points for a 
subsonic flow at M*.6 . The figure demonstrates the effect of 
the 3 successive iterations and also shows the convergence 
trend. The first iteration shows a considerable inaccuracy 



compared to the NASA code solution and the second and third 
iterations are in good agreement with the NASA solution. It is 
found that convergence occurs after 3 iterative calculations 
and additional iterations do not achieve better convergence. 
Thirty seconds of computer time on the CDC Dual Cyber were 
required to produce the result shown in Figure 4.1 . 




CHAPTER 5 
CONCLUSIONS 

Based on the cases examined, the computer programs 
performed well for two-dimensional airfoils of arbitrary 
thickness at moderate angle of attack. All the results were 
very accurate with the exeption of the Joukowski airfoil. The 
problem may be caused by the presence of the cusped trailing 
edge which introduces an additional condition. More likely, it 
is due to the error introduced in the conformal mapping. Hess 
[15] encountered similar difficulties when calculating the 
surface pressure distribution for airfoils with very thin 
trailing edges using a surface singularity distributions method. 
Hess used an additional parabolic vorticity variation that 
provided a satisfactory solution for thin trailing-edge airfoils. 
Using a similar method, Zedan [11] solved the direct and inverse 
problems of potential flow around an axisymmetric body using 
an axial source distribution. However, Zedan's method also had 
difficulty when solving the flow around airfoils with sharp 
corners or sudden changes in slope. 

In this study, the polynomial method provided an 
efficient and satisfactory solution to two-dimensional flow 
problems for airfoils with finite trailing edge angles. 

The use of double precision variables has proven to be 
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useful and gives stable and consistent results. Accurate 
solutions for the surface pressure distribution can be obtained 
on most airfoils by using 20 to 30 collocation points, 
especially if the calculation is made using the geometrically 
increasing grid. A typical case using 24 collocation points 
requires less than 1 0 seconds of computer time. One of the 
most important features of the method is its ability to deal 
with airfoils of any shapes by adjusting the value of the slope 
in the subroutine of the main computer program. 



APPENDIX A 

IDEAL FLOW OVER A JOUKOWSKI AIRFOIL 
UPWASH AND DOWNWASH VELOCITY CALCULATION 

A-1 Introduction 

The upwash and downwash velocities given by the 
equations 2.31 and 2.32 are obtained by approximating the real 
airfoil with an equivalent Joukowski airfoil of the same 
thickness. The solution for the flow about the Joukowski 
airfoil is accomplished using the traditional transformed plane 
and then the solution is shifted back into the real plane. The 
Joukowski method has the advantage that it determines the 
flow field anywhere in the real plane. In the present appendix 
we make an extensive study on the y-component of the velocity 
along the x-axis upstream and downstream of the airfoil as 
shown in figure A-1 . 

A-2 Flow about Joukowski airfoil 

Referring to figure A-2, we consider the transformation 
z=£+c 2 /£ from the transformed plane into the real plane, in 
which z=K+iy and £ = £+iii are complex variables. The 
transformation maps the circle of radius r Q centered at the 

origin of the t, plane into a Joukowski airfoil in the physical 
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plane, whose chord is sligtly greater than 2. In particular, the 
point £=1 is mapped into the sharp trailing edge of the airfoil. 
The shape of the airfoil is controlled by varying the two 
parameters m and 8 . For the present calculation, the airfoil 

becomes a symmetrical airfoil when 8= it. The radius r 0 of the 


circle and the angle p shown in figure 1 can be expressed in 
terms of m and 8. 

v 

u= 

These expressions will be used in the later analysis. The 

varialbies m and 8 which describe the displacement of the 

circle center in the £ plane are directly related to the camber 

ratio and to the thickness ratio by the following relations: 

— sin5 - 2 — 
c I 



(A.1) 
(A. 2) 


(A. 3) 


- — cos8 = -4=?- (A. 4) 

c 3^3 1 

where I is the total length of the chord of the airfoil. Finally 
we have a complete description of the airfoil parameters with 
c which describes the position where the circle in the £ plane 
cuts the £-axis and we note that c is given by: 
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c_ _ J_ 

I 4 

Under the same transformation, a uniform flow in the C, 
plane which makes an angle a with the horizontal x-axis, maps 
into a uniform flow with the same orientation in the physical 
plane. Let F be the complex potential of the flow in the £ 
plane, the complex potential consists of a uniform flow about a 
cylinder with the proper circulation that satisfies the Kutta 
condition in the real plane. In the current notation, the 
potential is: 

2 

F=U £ + Jl 

L c 

where: 

r - fp •§] ’ ia 
‘s l £ + me J e 

In order to satisfy the Kutta condition, the rear stagnation 
point needs to be positionned at an angle a+p which yields the 
relation: 

r 

sin(a+B) = — 

4nr U 
0 

The complex velocity of the flow about the airfoil is then 
derived by the equation: 


l 

+ i — ln-£- 
2 tt r 

0 


(A.5) 


(A. 6) 


(A.7) 


(A.8) 


C * 
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LV(z) * — 
dz 


the inverse of the Joukowski transformation is: 

C 




(A. 9) 


(A.1 0) 


We also have by definition Ul(z)*u-iu 

Inserting A.7 and A. 10 into equation A.6, F can be written as: 



After taking the derivative of F and separating the real and 
imaginary part, the y-component of the velocity along the 
x-axis upstream and downstream of the airfoil is given by: 



The mathematical problem for the equivalent Joukowski 
airfoil is now stated as follows: For a given airfoil with a 
specified thickness at a given angle of attack a , we have 
determined an approximation for the y-component of the 
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velocity upstream and downstream of the airfoil by using an 
equivalent Joukowski airfoil with same specified thickness and 
angle of attack. Note that the position of the leading edge of 
the described Joukowski airfoil is unknown. It may be 
calculated using the transformation and referring to figure 2, 
we have: 

C LE = b e m (A. 13) 

The parameter, b is given by: 

h.* 1 + inf cos(5-tt) -cos5] (A. 14) 

c c 

Using the definition of the transformation yields, 

in 

z=be +— (A.1 5) 


The x-coordinate of the leading edge is obtained by taking the 
real part of A.1 5. For example, for a symmetrical Joukowski 
airfoil with a specified thickness of 12% , we have: 

k L£ = -1.014405264 (A. 16) 



APPENDIX B 


AIRFOIL INTEGRALS 


The following are the Cauchy principal values for h z < 1 
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TT r n 1 1 (5) . . . (n-2) 
2 L 1_ J 2 (4) . . . (n-1 ) 
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16 . J „f „ 2 -l 

H-z, L 2 

17. [ i 3 V . 1-£ 2 ^ 

« 


-1 

1 


4 i 2 

H - — H 
2 


1 

8 


18. j d£ - -IT 

- 1 V 1 -$ <H-$) 
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+ * *************** ************** ★***♦***'************** ******'***'****** 

COMPUTATION OF POTENTIAL FLOW AROUNF TWO-DIMENSIONAL AIRFOILS 
USING A SINGULAR INTEGRAL METHOD: 

DIRECT PROBLEM 

** ****************************************************************** 
PROGRAM DBALCP (ANSWER .PLOT » OUTPUT * JOUKt. J0UK4* JOUK 1C* 
STAPE2=ANSWER ,TAPE6=J0UKC»TAPE9=JQUK4,TAPE1*=J0UK1C > 

** ** ************ * *** ******** ************ ************ * ** ************* 

DIMENSION X1C2 j5>*CPU(S35)*X2(235)*CP2C235»*X3C A13) »CP(41 0) 

S t PHUK2C5 1 *CC0M2 (2051 * C COM2 (2-1 5 )* CC0H3 ( 235 ) *CC0M4 ( 2C5 ) 

* ,PHU2f205>«GlSQ(2'5) 

DOUBLE PRECISION A<5 r*51)»OETER*ALPHA,Pl *X,ZfSX,R 


INPUT quantities: 

P0LYN0M OF DEGREE N 
N=2 0 


AIRFOIL: ItELLIPTlC AIRFOIL AT 12* ALPHA=D> 
2«NACAC0I2*ALPHA=:»IC0HP=! OR 2) 
3*NACACCi2,ALPHA=4,IC0MP=l> 
4*NACAC012 .ALPHA =1C*IC0MP=I> 
EyJOUKOUSKI AIRFOIL*ALPHA=C . JQUK? ) 

6, JOUKOWSKI AIRFOIL *ALPHA=4 * J0UK4 ) 

7t JOUKOWSKI AIRF0IL.ALPHA=1C. J0UK1S) 
8, CYLINDER > 


IMPLEMENTATION: 

IP0S=2 


i .E qually spaced points 
* .EXPONENTIAL QUADRATURE 


IPCLY=1 


polynqh: 


lyPOUER SERIES OF X 

2y SQUARE ROOT AND POWER SERIES OF X 
3*AIRF0IL POLYNOMIALS 


ICQMP=1 INCOMPRESSIBLE 


IC0MP=2 


::compressible 


EXPONENTIAL RATIO 
R=1.1D+C 


ANGLE OF ATTACK 
PI=3 .141592654D*CG 
ALPHA=C- «D* £ 

IF(KAI.£Q.3 .OR. KAI.EG.6) 
IF (KAI.EQ.4 .OR. KAI.EW.7) 
KA2 = 1 


ALPHA =( 4.D*0)*PI/(i8C.D*0) 
ALPHA = ( 1 C.D*0)*P 1/(18 C.D*u ) 


COMPUTES THE MATRIX ELEMENTS ( N EQUATIONS, N UNKNOWNS) 


ORIGINAL PAGE IS 

OF POOR OTJatt'tv 


non non non oonoo onnno 


continue: 

1=1 

CALL POSIT (X.I.NtIPOStRJ 
CALL MAT1C A*I*X*IPOLY,N) 

CALL DFDXCZtKAI tXt KA2) 
IFCDABSCZ>.LE.C1.D-5)J THEN 
N=N-1 
GO TO 33 
END IF 

CALL MAT2(A*Ntl t XtIP0LY»KA2tZ) 


COMPUTES THE UPWASH AND DO«N«ASH VELOCITY INTEGRALS 

sx=e*o*o 

IFtALPHA.EQ. CO. □♦'>>> 6C TO 27 
CALL VELEXT (SX«ALPHA»X> 

27 CONTINUE 

IFCKA2.EG.1) AC I« N+1)=-SX-QC0S (ALPHA )*DS IN C ALPHA >* 

* DLQGCCCi.O*OJ*XJ/C (l.u«-0>-X>>/PI 

IFCKA2.EQ.2) AC It N+l)=-SX-*-DCOS C ALPHA! ♦OSIN C ALPHA )* 

S DLOGCCCl.D+0>*X>/( C l.D«-0>-X))/PI 

I=I^1 

IFCI.LE.N) GO TO 2 " 


GAUSS-UORDAN REDUCTION 

DETER=C!.D-C 

CALL GAUSS C A tN. DETER ) 

VRITE C2t 2T 3) of TER 
URITEC2t2i8> 

URITEC2t-201) C A (It N+l) * 1=1 tN) 


COMPUTES THE PRESSURE DISTRIBUTION USIN6 THE MATRIX COEFFICIENTS 
PHU1=X COMPONENT OF THE DISTURBANCE VELOCITY 


XI C 1 >=~i •♦2./FLOAT (201) 

DO 30 I=lt 2 CC 
X=DBLECX1 Cl) ) 

CALL DFOXCZtKAItXfKAZ) 'i 

IF CDABS C Z ).LE«( 1. D~E > > THEN 
CPUCD=CPUCI-i> 

GO TO 36 
ENDIF 

CALL CPOICAtZtNtPHItXtIPOLYJ 
PHUl(I)=PHI-l. 

CPU(I)=1.-CPHI**2>*C1.*SNGLCZ)*SNGLCZ>> 

Q1SQ ( I) =(PHI**2 )*( 1 .♦SNGLC Z J*SNGLCZ)> 

16 Xl(I+l)=XiCI>*2./FLCAT(2Sl> 

3~ CONTINUE 

CALCULATES THE COMPRESSIBILITY FACTOR 

IFCIC0MP.EG.2) CALL CALCOMC XltPHUl t PHU2 tCC0HltCC0M2 tCC0M3 1 
$ CCQK4 tKA It KA2) 

«RlTEC2t22:) 


ORIGINAL PAGE IS 
OE POOR QUALITY 
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DO 40 I=2»2C0*5 • 

4S MRITEC£*.221> XI (I) *PHU2 ( I> *Q1SG(I) » CCOHi (I ) *PHU2 C I > 

220 FQRMAT(/7X,*X'*13X* *PH1 • *5X * *G1 SQUARE* » 3X* *CPI C0HP**7X, 

* 'PHI'/) 

221 FORMAT (5 (Fll *5, 3X) > 

IF( ALPHA *£Q a O.D+C ) } GO TO 31 
IFCKA2.EQ.Z) GO TO 31 


FOR A NON ZERO ANGLE OF ATTACK* CALCULATES THE 
LONER PLANE DISTRIBUTION 

KA2=2 

DO 32 1=1, 2CC 
X2(I)=X1( I) 

32 CP2C I)=CPUCI > 

GO TO 33 
31 CONTINUE 


READS NUMERICAL PRESSURE DISTRIBUTION FOR COMPARISON 

CALL C0MP(N2,CP,X3*KAI,IC0MP) 

223 FORMATl/, 'DETER =',Cl*.12/> 

2 21 F0RMAT(lH,13Di; .7> 

2 28 FORMAT t/ 'POLYNOMIAL COEFFICIENTS'/} 

2 25 F0RMATC2F1^.6) 

2 26 FORMAT (F14«6) 


PLOTS THE PRESSURE DISTRIBUTION 
IF(ICOMP.EQ.l) THEN 

CALL FI PLOT l XI, CPU ,X2*CP2 , X3 ,CP *N2 • KAI, ALPHA *IP0S* IP0LY*N,R > 

CALL COMPLCT (XI *CCCK1*CC0MZ *CC0M3* CC0M4 * X3 *CP,N2 *KAI* ALPHA* 

* IPOS* IPOLY * N*R ) 

ENOIF 

STOP 

END 


**« *♦***•*■***■*■•■•*** + ***■***♦**** »★■**■* ***« * ***"*#*4r * ***** * *** ** #*** 

SUBROUTINE F IPLOTCX1 ,CPU*X£ ,CP2,X3,CP,N2 * KAI * ALPHA, IPOS , IPOLY, 

* N,F. } 


PLOTTING SUBPROGRAM 

DIMENSION Xl C21 5), CPUC2 05 )*X2(2D5)* CP2C2S5)«X3 (4 13 >*CP(4l") 
DOUBLE PRECISION ALPHA, R 


CALL PLOTS CC,0,4LP LOTI 
CALL FACTOR C } 


SCALES OF THE AXIS 
Nl=2"0 

IF < IPQLY»Eu«3) Nl=194 


ORIGINAL PAGE IS 
OE POOR QUALITY 
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IFCKAI.EQ.l .OR* KAI.EQ.8) N1=1QS 
XiCNl*i)=-l. 

X1CN1^2) = .4 

IFCKAI.EQ.l .OR. KAI.EQ.8) XlCNl+2)=.2 
X2CN1*1)=X1CN1*1> 

X2 C N1+2)=X1 C Nl+2) 

X3 CN2+1) =X1C Nl+1) 

X3CN2«-2)=XlCNl-*2) 

CPUCN1+1)=1. 

CPUCN1«-2)=CPCN2*2) 

CP2 CNld )=1. 

CP2tNl*2)=CPCN2-»‘2) 

CPlN2*l)=l. 



AXIS AND ORIGIN 

CALL AXXSC 2.2*2 .2* 1HX*-1 * 5 • #u.l #X1C N1-*-! ) *X1 CN1+2 ) > 

CALL AXISC2.2.£.2f2HCP*2»G.C»9C.3*CPUCNl^U*CPUCNl-»-2)> 
CALL AXISC2.2*8.2*2H #-l *5. *0.*999 • *1.) 

CALL AXISC7 .2*2 .2* 2»£*6.*93«*999.»1 •) 

CALL 0RIGINC2«2*2.2» £) 


DESIGNATION OF THE CURVES 
CALL PLOT C .5*5.75*3) 

CALL PLOTC .8*5.75*2) 

IFCKAI.EQ.l) CALL SYMB0LC1 . #5.75 *. I *»EXACT ANALYTICAL *# 0. *1 
IFCKAI.EQ.f .OR. KAI.EQ.6 .OR. KAI.EQ.7) THEN 
CALL SYMBOLC I .#5. 75# • 1# a THEORETICAL* • 0.# 11) 


ENDIF 

IF CKAI.EQ.S) CALL SYMBOLC 1 • #5.75 ».l » ^THEORETICAL • • 
IF CKAI.EQ.2 .OR. KAI.EC.3 ) THEN 

CALL SYMBCLtl.»5.75*.l# • NUMERICAL TAB L123**C.* 


3 . # 11 ) 

18) 


IF CKAI.EQ.4) CALL SYHBOLC1 • #5.75 #.l # 'NUMERICAL DATAL133** 
Sl.tlS) 


6) 


CALL LINECX3#CP»N2#1#:#2) 

CALL LINECX1#CPU#N1#1#-6»1) 

IF C ALPH A.NE .CD. D+C C ) ) CALL LINECX2*CP2*Nl*i*-6*l) 


CALL SYMeOLC.6»5.5».l*I*C.*-l) 

CALL SYMBOLd^.^#.!#* POLYNOMIAL N=»*3.#13) 

CALL NUMBER C2.5»5.5*.1#FL0ATCN)#G.#3) 

IFCIP0S.EQ.2) THEN 

CALL SYMBGLC4.#5.75*.1* , R= , »C.#2) ' 

CALL NUMBER C 4.3 #5 .75# «1#SNGL(R) • 0.»1) 

ENDIF 

IFC IP0LY.EQ.2) CALL SYHBCLC4.*5.5#.1#»POLY=2*»0. *6) 
IF C IPOLY .EG. 3) CALL SYMBOL C4. #5 .5* .l.^POLYsS*#^. #6) 


IFCKAI.EQ.l) CALL SYMBOLC. 2*6. 25#. 14#»ELLIPSE ALPHA = 3 • * 0. * 16) 
IFCKAI.EQ.fi) CALL SYMBOL C .2 #6.25 #.14# •CYLINDER* • 3. *8) 
IFCKAI.EQ.2) CALL SYMBOLC. 2*6. 25#.14 
S»* NACA 3 012 AIRFOIL ALPHA=0»»0. #25) 

IF CKAI.EQ.3 ) CALL S YM60L C . £ *6.25 #.l 4 
*»• NACA 0012 AIRFOIL ALPHA=4 *# 3. # 25 ) 

IFCKAI.EQ.4) CALL SYMBOLC. 2*6. 25*. 14 
Sa'NACA G012 AIRFOIL ALPHA=1 5 •# ?• *2 6) 

IFCKAI.EQ.5) CALL SYMBOLC. 2 *6.25 #.14 
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St* JOUKOUSKI AIRFOIL ALPHA=0* , 0 . • 25 1 
IF C KAI *£Q« £ 1 CALL SYMBOL C .2 *6. ?5 # . 14 
S« 'JOUKOUSKI AIRFOIL ALFHA=4 * « S . » 25 1 
IFCKAI.EQ.7) CALL SYMBOL C.Z to. 25 *. 14 
S.* JOUKOUSKI AIRFOIL ALPHA=10**S. t2fc) 
CALL PL0TC-,5t-.5t999) 

7 CONTINUE 
RETURN 
END 


ORIGINAL' PAGE IS 
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** ********************* ************************************ ********* 
SUBROUTINE DFDX CZ,KAI,X, KA2 ) 

** ********* ********************************* ****************** ****** 
SUBPROGRAM THAT CALCULATES THE DERIVATIVE OF THE AIRFOIL 
DOUBLE PRECISION Z«X 


ELLIPTIC AIRFOIL 

IFCKAI.EQ.il Z=-C.12DfC>*X/DSaRTtlI.D«-ei-X*X) 

CYLINDER 

IF (KAI.EQ.&l Z=-X/DSGRTCtl.D^il-X*X) 

NACA 1?12 

. IF CKAI.EQ.2 .OR. KAI.EG.3 .OR. KAI.EQ.4) THEN 
Z=(1.2D-*C)*t (.1 74225 C 1 *”! l/DSQRTC iX*-i 1.C+0 11/ 

* C2.D«-0)l-C.e6ZD*C>-C.35i£DO>*C CX*<1.D*Q>)/<2.D«-C>) 

S ♦C.42645D-*-01*C CtX«-Ci.D*-3>>/C2.D*0)>**<2.D-*-G > > -C .2 73D+C >* 

S CC 4X*Ci.D+0>>/C2.C*C)l**C3.D^C>l) 

end if 

JOUKOUSKI AIRFOIL 

IFCKAI.EQ.5 .OR. KAI.EQ.6 .OR. KAI.EQ.71 THEN 
Z=C C4.D-01/C (3.D+C )*DSCRTC3.C*">)>*C.12D*01* 

$ C-DSQRTlCl.D*01-X*X>-tCl.D+‘ 1 -X 1 *X/DSQRT t Ci.D+O >-X*X ) ) 

END IF 

IF (KA2.EQ.2) Z=-Z 

RETURN 

END 

** *************** ****************************************** ********* 
SUBROUTINE VELEXTCSX *AL2*X) 

***************** ****************************************** ********* 

CALCULATES THE UPUASH AND DOUNUASH VELOCITIES 

DOUBLE PRECISION VEL *XLt AL2 *RM « XSLE # SX*X*PI 
PI=3.14* 5926 s 4S + *j* 

RM= COEFFICIENT FOR THE UPUASH AND DOUNUASH JOUKOUSKY 

CALCULATION OF THE EXTERNAL VELOCITY 
RM=SQRTCTHICKNESS*THICKNESS/27*C AMBER *CAHBER/4) 

RM=. 0461880 2D* 1 

XSLE=-1 • 0144 052640*7 

XL=-C1.D*0>-C.02D*C)/C2.D*0) 


UPUASH CALCULATION 
DO 25 IL=i»5 00 
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XS=C1.D*D>-C tl.D+C>-XL>*CC2.D*5)-XSLE)/t2.D+C> 
VEL=CDSfiRTC CXS-Cl.D*G))/(XS«>Ci.D-*-G>>)*< ( txs/C2.D*Q>- 
SDSQRT CXS*XS/ (A.D*-!)-* .250+C )>)/ CXS/ C2 .D* 3> -QS3RT CXS*XS/ 
$(4.D+0)*(.25D+0))*RH))**(4.D+!2>)-{ 1.0*P> )*DSIN< AL2> 
SX=SX-»-VEL*C.32DfS)/tX-XL)/PI 
25 XL=XL-t.22Q«-3> 


DOUNU ASH CALCULATION 
XL=tl.D«-C>-*t .D2D*S)/(2.0^C1 

DO 26 IL=1 *5 DO 

XS = < 1 •□♦3) -C C1.D+0>-XL>*«C1.D^1)-XSLE)/C2.D*0) 
VEL=CDSeRTCCXS-CI.D*0>J/CX£*Cl.O*:J))*C<CXS/C2.D*D>**‘ 

SDSQRT CXS*XS/CA •0'*-G »-< «25D^0 > > J / l XS/ <2 *D^0> ♦0S3RT CXS*XS/ 
tCA.D+Q)-C.25D^: )>*RM)>**l2.D^3J)-tl.O-*-C) )*DSINCAL2> 
SX=SX«-VEL*C. 02D*-3)/lX-XL)/PI 
26 XL=XL-*-{. C2Q'*' G) 

RETURN 

END 

*■* ,**•*♦■****■.**♦■***-***■****•.♦******♦**■*■.*.•■*■*************•* ************ 
SUBROUTINE 6 AUSS CA ,N »CETER > 

** ***************************************************************** 

GAUSS-JORDAN REDUCTION 

THIS SUBPROGRAM FINDS THE SOLUTION VECTOR CORRESPONDING TO A 
SET OF N SIMULTANEOUS LINEAR EQUATIONS USING THE GAUSS-JORDAN 
REDUCTION ALGORITHM WITH THE DIAGONAL PIVOT STRATEGY. 


DOUBLE PRECISION A 15 D»51) 1 DETER *EPS 
EPS=1.D-1J 
NPLUSM=N-M 

C BEGIN ELIMINATION PROCEDURE 

DETER=U.O*3| 

DO 9 K=1,N 

C ..... UPDATE THE DETERMINANT VALUE •••• 

DETER=DETER*ACK*K) 

C ••••• CHECK FOR PIVOT ELEMENT TOO SMALL 

IF(CABS(A(K 9 K>>.GT.EPS) GO TO 5 
GO TO 7 

C ..... NORMALIZE THE PIVOT ROy ..... 


5 KP1=K*1 

DO 6 J=KP1,NPLUSM 

6 A( K , J)=A (K , J )/ A (K, K) 

ACK,K)=i.D*C 

C ELIMINATE KCTH) COLUMN ELEMENTS EXECPT FOR PIVOT 

00 II 1=1, N 

IF C I.EQ.K .OR. A*I,K).EG.i O.D-*S> ) SO TO 11 
DO 8 J=KP1*NPLUSH 
3 ACI, J)=A(I,J)-A(I,K)*A(K,J> 

AtI,K)=C.D*0 
11 CONTINUE 
9 CONTINUE 
GO TO 1 

7 CONTINUE 
URITEC2»2D2) 

2C2 FORMAT*/, 'SMALL PIVOT - MATRIX MAY BE SINGULAR*) 

STOP 

1 CONTINUE 
RETURN 
END 
C 
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C ** ****** ♦********+♦****'**+******'** ****'•''» ***-***** ******* **** A"******* 

SUBROUTINE POSITlX,I # N,IP0S,R> 

** ************* ******** ************* **** *************************** 

SUBPROGRAM THAT IMPLEMENTS THE LOCATION OF THE SET OF QUADRATURE 
POINTS 


LIST OF THE VARIABLES: 

XP LOCATION OF THE POINTS 

N NUMBER OF POINTS 

IPOS TYPE OF SPACING 

R EXPONENTIAL RATIO 

DOUBLE PRECISION XPC 5 OJ * R*DELTA tDX ,XL.X 
IFUPOS.EQ.i J THEN 

XPtlJ=-tl.D*00J«-t2.[j*0e J/DBLEtFLOATtN+U J 
00 10 J=2 ?N 

XPtJJ=XPlJ-iJf€2.D*:CJ/DBLEtFL0ATtN-*-lJJ 
1' CONTINUE 
END IF 




IFtIP0S.EQ.2J THEN 
NTE=N/2 
NLE=N-NTE 

D ELTA=t l.D+C J/( C C 1.D+ 0 J- tR**NLE J J/ttl.D+OJ-RJ 


XL=-1.D+ 


DO 1 J=1 


*NLE 


DX=DELTA* CR**( J-l J) 


XPtJJ=XL«-LX/t2.C-*-0J 

XL=XL+DX 


> 


DELTA=tl.O-K )/t ttl.D«-CJ-tR**NTEJ J/ttl.O«-0J-RJ J 

xl=i.c+: 

DO 2 J=1*NTZ 
DX=DELTA*tR**tJ-l JJ 
XP(N-*-l-J)=XL-DX/C2.C*0) 

2 XL=XL-OX 
END IF 
X=XPtI> 

RETURN 

END 


SUBROUTINE MAT1 CA«IsX*IPOLY *NJ 


DOUBLE PRECISION A (5 1*51 ) ,X,RMT ,PI ,RM1 
PI=3 .14 15926 540 •♦■Q 

SUBPROGRAM THAT COMPUTES THE MATRIX ELEMENTS 
IFCIP0LY.ES.1J THEN 

At I*U = tDLOGttt l.D*: J*XJ/t Cl.Df OI-XJ JJ/PI 
DO 1 J-2 y N 

AtI*JJ=X*ACI»J-lJ“tt 1 •D‘*’C 0 J—DBLE CFLOATt l**l J** IJ-1J J JJ/ 
* tDBLEtFLOAT td-i J J*PI J 
1 CONTINUE 
END IF 

IFtIP0LY.EQ.2J THEN 
A 1 1 * 1 J=X 


c 

c 

c 


ACIt2)=X*X-«.5D+0) 

A 1 1 1 3)=X* At 1*2) 

ACIt4)=X*ACIt3>-Cl.D«-C>/t8.Q*0> 

AtIt£)=X*AlI*4) 

D0 


original page is 

OR POOR QUALrry 


DO 3 IMT=i» J-3 »2 

RHT=RMT*DBL£CFLOATCIMT/tIMT-*'l )) i 
A( I,J)=X*A(I* J-l)-«.5D*-C)*DBLEtFL0ATCl-t t-l>** t J-l) ))) 
* *RMl/D8LE CFLOAT < J>> 

continue: 

END IF 


IF ( I POLY *EQ • 3) THEN 

Ati,i>=i.c*-: 

AtIt2>=Cl.Df2>*X 

AClt3)=(.5D+G)*X^X*X 

A(I*4) = t.5D*D)*t.5C-*-0>*X*X*X*CX**2) 

DO 4°J=E,N 
RMT=1.D*? 

RKl=l.D«-v 

DO 5 INT=l,J-3 ,2 

5 RHT=RMT*DELECFLGATCIMT/CIMT*i>>> 

DO 6 IHT=1 » J-2 *2 

6 RM1=RM1*D6LECFL0AT(IHT/CIMH»1>>> 

A( I. J)=X*A(I, J-i)^C.5D^C)*DBL£CFL0ATCl-CC-i) **IJ-1> )) ) 
$ *RHT*« .5 d->2>*DBL£CFL0AT<I-CC-1>**CJ)>))*RH1 

4 CONTINUE 

ENDIF 


RETURN 

END 


C 

r 


** *************************** * ******* *** * ******************** ****** 

SUBROUTINE MATE (A»N»I»X»IP0LY»KA2»Z) 


DOueLE PRECISION A(5 :*£!>» X*Z*RS 
C 

RS=1.D*0 

IF IKA2.EG.2) RS=-1*G+C 
C 

IFCIP0LY.EC.1J THEN 
A(I*l>=ACI*i)-RS/Z 
DO 1 J=2tN 

AtI*J>=AtI*J>-RS*(X**t J-UJ/Z 1 

1 CONTINUE 
ENDIF 

C 

IF( IPOLY.EQ*Z) THEN 

A(I»1>=A< I*l)-RS*DStRTC ( 1*D*3 )-X*X>/Z 
DO 2 J=2#N 

A< I*d)=AC I» JJ-RS*DSURTt C1.D^3>-X*X>*C X**U-1 ) )/ Z 

2 CONTINUE 
ENDIF 

C 

IF ( IPOLYaED* 31 THEN . % , 

ACItl)=AC I*i)*RS*DSORTC< Cl.D*C>*X)/( « 1*D*3)-X> )/Z 
DO 3 J— 2 t N 

ACI*d)=ACI*d)*RS*DSaRTCtll*D^3)^X)/t Cl.D«-0>-X> > * IX**C J-l > >/Z 

3 CONTINUE 
ENDIF 


on o non on on o non on 


R r TURN ORIGINAL PAGE IS 

END OE POOR QUALITY 

** ********************* ************************************ ******* 
SUBROUTINE CPDI (A,Z,N,PHI,X,IPCLY> 

* * ************* *************************** ** ********************** 

PRESSURE DISTRIBUTION 

DOUBLE PRECISION A(5C,51),Z,X 

IFCIP0LY.EG.1J THEN 

PHI=SNGHAC1,N*I) )/SNGL (Z) 

DO 1 J=2,N 

„ PHI=PHI*SNGL{ A C*J* N*l) )*SNGL( X** ( J— 11 l/SMGL (Z 1 

1 CONTINUE 
END IF 

IFCIP0LY.EQ.21 THEN 

. PHI=SGRT ( l.-SNGL( X*X 1 >*SNGL( A(l, N*i) l/SNGLCZ) 

DO 2 J— 2»N 

PHI=PHI*SGRT(1«-SNGL(X*X))*SNGL( A(J,N+U 1*SNGL(X** C J-l> )/ 

S SNGL(Z) 

2 CONTINUE 
END IF 

IFCIP0LY.EG.31 then 

PHI=-SGRT ( (1 .+SNGLCX ))/(l.-SN6L(X) >) *SNGL ( A( l»N*iJ)/SNGL(Z) 
DO 3 d=2»N 

PHI=pHI-SQRTC(l.-*SNGL(Xl l/(l.-SNGL(Xlll*SNSL(A(J,N+ill* 

S SNGL (X** ( J-l 1 1/SNGL(Z 1 

3 CONTINUE 
END IF 

RETURN 

END 


SUBROUTINE COMP IN2.CP.X3 .KAI.ICOMP) 


READ THE CP DISTRIBUTION FCR COMPARISON 

DIMENSION CP (41 01, X3 (4101. XP1C5 01. CPI (5 0 
$ ,XP2(401,CP2(4C 1,XP4(401,CP4(4S1 

* ,XP3C 701,CP3(7C1 


DATA 

FOR NACA .012 

ALPHA=0 



DATA 

XP1/1.CCC0 00* 

•997865, 

.991852, 

•982340, .969515, 

$ 

•953555. 

•934629, 

.912919, 

•888573, .861836, 

% 

.832804. 

•802769, 

•768917, 

•734470, .698658, 

* 

.661720, 

•623698, 

.535439, 

.546589, -537593, 

S 

•468694 , 

• 43.' 128, 

•392122, 

•354896, .318658, 

% 

.283605, 

.245925, 

.217791, 

•187366, .158834, 

% 

• 132244, 

•1C 7820, 

•385651, 

.165848, .248511, 

$ 

$ 

•133727, 

• C 00100/ 

• w2i5 71 , 

.212098, 

•3Q5346, .001327, 

DATA 

cpi/i.cco:co. 

• 35-/876, 

•26252 : , 

.203777, .158301, 

$ 

.119665, 

• 166223 , 

• 05 6265 , 

•128921, .003559, 

% 

-.020291, 

-.143022, 

-•<^64972 ,— 

• 286434 ,-.l 07666, 

S 

-.128882, 

-.15:253, 

-.171887,- 

.193848, -.216129, 

s 

-.238655, 

-.261284, 

-.283815 

.305942, -.327360, 


********* *» «*«* «*«*<*«* «*«*«*«***** *♦*.***,***» ******** 



DATA FOR NACA '01 
DATA XP3/1.C33C0C 
• .929474 

•751896 
•517355 
•283614 
. 1 : 2:67 
.008384 
.018953 
.132250 
.327623 
.566061 
.793736 
.953575 


ALPHA=4 
.996731* 
.901:76* 
.7 07743* 
•468706* 
.241747* 

• 0 75*t52 * 
•C 02081* 

• 033731 * 
.165765* 
.3 7340 8, 
•614345* 

• 832222 * 
.973:43* 


98 754 0*. 97 3243*. 95357 5* 
863734* .832822*. 793736* 

661736. . 614348* .56606 1, 
42Q57 6*. 373*0 8*. 327623, 
2D2362*. 165768* .132250, 
C52614*. 033731* .01395 3, 
003:33*. 002081* • 008384* 
C52614*. 075452*. 102067, 
202262*. 241747*. 283614, 
423576*. 468736*. 51735 5* 
66173 6*. 7 3774 3*. 751896* 
868734*. 9D1376*. 929474, 

987543. . 996731/ 


DATA CF3/1.D64 172 
. 13*079 
. 212:01 
-.063249 
-.083448 
.126406 
.981974 
-1.813338 
-1.1341*66 
-.685359 
-.364722 
-.124925 
.118939 


•374255* .280179* .218254, .171865* 

.1 22582* .075460, .251872, .030885, 

-•C 2 5333 » -.021*03, -.036427, -.050363* 
-.:740C6* -.082602* -.087861* -.388644, 
-.17C629, -.048107, -.313241, .036378, 

•21227, .329655, .523111, .730817, 

1.021316, .335786* -.951149,-1.685110* 

-1. 693127, -1.531 778,-1. 38 1722, -1.25 3 222 * 
-1 .229778 , -.934240* -.84572? * -.762944, 
-.612491, -.544249* -.480336, -.420631* 
-.312287, -.262711, -.215467, -.169755, 
-.:eOC4G* -.334382* .313221* .063638, 

.181296* .257856* .36468CV 


DATA FOR NACA 012 ALPHA=1 " 

DATA XP4/-. 9944,-. 989C,-. 9833, -.9690, -.9444, 
-.9033*— .SCCD * — .7000*-. 60 :0*-.4?3G, 

— .2301* • v 0 GO * .2003* .43-0* .63*20* 
.8003, .8611, .9778* .8611* .8333* 

.6 ?•::*, . 4001 * .2000* . OC 0* — .2 .?■? - * 

- • 4 0 0 0 * - • 60 C C* » - • 8 2 G 0 * - • 85 C 0 * -• 9 44 1 * 

-.9 722, -.9801* -.9910/ 

DATA CP4/-6 .11 . 0,-6. 4166 , -6.0C 30 ,-5.0000,-4.3000, 
* -3.303* 0,-2. 1670, -1.7222 *-1.4400,-1. 322 0 » 


c 

c 


c 

c 


c 

c 


c 


c 

c 


c 

c 


% 

s 

s 

s 

s 


- • 75 0 C , 

-• 055 5, 
•1550* 

• 3 3 5 C * 

• sr-'c* 


.5670, —. 4440 , 

• 5Sro t .2778, 

•1667, .1667, 

• 4170* .687 Z * 

• oe.c,-i.co: •>/ 


-.3555, 
• 1 667 * 
. 2220 , 


• 194Q 
1667 


.25* 0 , 


.7780, 1*0020, 


(AIRFOIL NACA12 AT 4 DEGREES OF ANGLE OF ATTACK) 

IF (KAI.EQ.3) THEN 

N2=64 

CPC 66 )=- *5 
DO 41 1=1,64 
X3 ( I ) =XP3( I) 

CP( I )=CP3( I) 

CP(I)=CPCI)*SQRT(i.-.5*.5> 

41 X3CI)=2.*(X3CI>-.5) 

END IF 


.ORIGINAL PAGE IS 
OR POOR QUAI.ITY 


(AIRFOIL NACA 0 C 12 AT II DEGREES OF ANGLE OF ATTACK) 

IF (KAI.EQ. 4 ) THEN 
N 2=33 

CP ( 35 )=-l • 

00 48 1 = 1,33 
X 3 ( I)=XP 4 ( I ) 

CP ( I ) =CP 4 ( I ) 

48 CONTINUE 
END IF 

(AIRFOIL NAC All AT C DEGREE OF ANGLE OF ATTACK) 

IF (KAI «EQ .2 .AND. ICGHP.EQ.l) THEN 
N2=41 

CP (N 2+2 )=-• 3 
DO 44 1 = 1 1 N 2 
X 3 ( I ) =XP 1 ( I ) 

CP(I)=CP 1 ( I) 

44 X 3 ( I ) = 2 .*( X 3 (I) -. 5 ) 

END IF 

IFCKAI.EQ .2 .AND. ICOHP.EQ.Z) THEN 
N 2=33 

CP ( N 2-*>2 ) =- • 2 
00 45 1 = 1 , NZ 
X 3 ( I)=XP 2 ( I) 

CP ( I )=CP 2 ( I ) 

45 X 3 ( I )=£•*( X 3 ( I) -. 5 ) 

ENDIF 

(ELLIPTIC AIRFOIL) 

IF (KAI.EG.l) THEN 

N2=1C1 

CP ( 10 3 )=-« 2 

X 3 (l)=- 1 . 

CP( 1 )= 1 . 

DO 46 1 = 2, 101 
X 3 CI)=X 3 (I-l )«- 2 ./irC. 

Q£ = (l.«-.I 2 )/SGRT(l.-*-. 12 *. 12 *X 3 ( I)*X 3 (I)/( 1 .-X 3 (I)*X 3 (I ))) 

46 CP(I)= 1 .-QE*GE 
ENDIF 


(CYLINDER) 

IF (KAI.EQ.3) THEN 

N2=101 

CP (11 3)=-* 8 


I 


no on no no nnoo no n on no on 


CP(1 )=1. 

DO 43 1= 2*101 

X3( I ) =X3CI-i >♦!•/! *C. 

43 CP(I)=-3.^4.*X3CI)*X34I> 

ENDIF 

(JOUKOWSKI AIRFOIL AT I DECREE! OF ANGLE OF ATTACK) 

IF(KAI.E0*3) THEN 
N2=2'22 

CP l N2+2 ) =- • 3 
DO GO 1=1, 2C2 
63 READ <8*2 3 d) X3CI)tCPCI) 

END IF 

(JOUKOUSKI AIRFOIL AT 4 DEGREES OF ANGLE OF ATTACK) 

IF(KAI*EQ*6> THEN 
N2=434 

CP(N2«-2)=-*5 
DO 61 1=1,202 

61 READ (9,205) X3(I),CP(I) 

DO 62 1=1*232 

62 READ (9*235) X3(4C5-I ),CPC405-I) 

END IF . 

(JOUKOySKI AIRFOIL AT 13 DEGREES OF ANGLE OF ATTACK) 
IFCKAI.EQ.7) 

N2=4 04 

CP(N£+2)=-l. 

DO 63 1=1, 2C2 

63 READ ( 10*235) X3(I),CP(I) 

DO 64 1=1,202 

64 REAC ( 10,.203) X3 (425-I)»CPC4C5-I> 

END IF 

2 C5 FORMAT (2F1C. 6) 

RETURN 

END 

** a************************-************************** * ****♦■»** * 

SUBROUTINE CALCOMC XI ,PHU1 * PHU2* CCOK 1 * CC0H2 * CC0M3 ,CC0M4 ,K AI , KA2) 
** ************************************************************* 

CALCULATION OF THE COMPRESSIBILITY FACTOR 

DIMENSION XI (2:5)*PHUK205)*PHU2(2( 5) ,PHU3 (205) ,PHU4(2 35) 
DIMENSION CCGMl (2C5) *CC0H2 (2 35) * CC0M3C2C5) *CCQM4 (2 25) 

DOUBLE PRECISION A(5 C*51 ), DETER ,Z 
N1=2CG 

PI=3. 141592654 

PHU1 =X-COMPONENT OF THE DISTURBANCE VELOCITY 
CC0M1=PRESSURE COEFFICIENT 


CALCULATES CP(IST) 

CALL CALCPC (XI, PHU1 * CC0M1) 

CALCULATES CP(2ND) 

DO 1 1=1,2 32 
PHU2 ( I ) =PHU 1 (I) 

C CALCULATES THE POLYNOMIAL THAT FITS BETWEEN -1 AND ♦ 1 


non non no n no n no 


CALL MATCQMCA*X1*PHU1J 
CALL GAUSS (A *2C *DETER) 
RINT=ALQGttl.«-XltI>>/<i.-XlCI>l>/PI 
PHU2 CI>=PHU2 CD+SNGLCAC 1*21)>*RINT 
DO 2 2M=2§ 2 0 

RINT=RINT'X1 tl>~Cl.-FLDATlC-l>**CIM-l>)> 
t /(FL0ATCIK-2)*PI> 

2 PHU2CI>=PHU2 Cl) ♦SNGL CACIK*2D) *RINT 

i continue: 

CALL CALCPCCXltPHU£«CCCM2) 


0*lGlH Ar n 


CALCULATES CPC3RDJ 
DO 3 1 = 1 * 2 :? 

PHU3 < I)=PHU1 Cl ) 

CALCULATES THE POLYNOMIAL 


THAT FITS BETWEEN -1 AND 


CALL HATCQKt A*X1 »PHU2) 
GAUSS (A *20 *DETERJ 


CALL 

RINT=AL0GC U.*XlCI)>/t2.-XiCI>)>/PI 
PHU3C I)=PHU3 CI> ♦SNGLCAC 1* 21 )) *R INT 
DO 4 IH=2,20 

RINT =RINT*X1 Cl >“C 1 »“FL0AT C C-l )** C IM-ll ) 1 
J / 1 FLO AT CIM“1> *P I) 

4 PHU3 tI)=PHU3CI>+SMGLtACIM»211>*RINT 
3 CONTINUE 

CALL CALCPCC X1*PHU3 *CCCM3) 

GC TO 7 


CALCULATES CPC4TH) 

DO 5 1=1 .2 TO 

PHU4 ( I J =PHU1 Cl) , „„ . . , 

calculates the polynomial that fits between -a and *1 

CALL MAT COM < A*X 1«PHU3) 

CALL GAUSSCA*2' *OETER> 

RINT=ALQGC C1*+X1C1))/C1 »-X 2 Cl) ) ) /P I 
PHU4 C I)=PHU4 Cl ) +SNGL CAC 1*22)) *R INT 
DO 6 IH=2,2? 

RINT =RINT*X1 C I >-C 1 « -FLOAT C C-l)** CIH-l))) 

PHU4(I)=PHU4CI)+SNGLCACIH*2D)*RINT 

CONTINUE 

CALL CALCPCC XI »PHU4*CC2M4 ) 

CONTINUE 

RETURN 

END 


** •*<* *** Ifr***-* *****★♦*♦♦ * *** ********** ***»************?* ******* 4 


SUBROUTINE CALCPCCX1 *PHU1 *CC0H1 ) 


CALCULATES CP FROM PHI 

DIMENSION XI C2?5)*PHU1C205) »CC0M1C2Q5> 
DOUBLE PRECISION Z*X 


PHU1=X-C0MP0NENT OF THE DISTURBANCE VELOCITY 


Ni=2C2 
KA 1=2 
KA£ = I 

DO 1 1=1 *2 CC 

X=DEL£C Xi Cl ) ) 

CALL DFOX CZtKAI *X »KA2) 
IFtDABSCZ).LE.Cl.D-5)> 
CCQHl C I)=CCOKi CI-D 


THEN 


noon on on nonn no 


GO TO 2 
ENDIF 

Q2=CC1.«-PHU1CI)>**2>*C1.«-SNGLCZ>*SNGLCZ>> 
RGAM=1.4/Ci*4-1.) „ 

RFAC=1. +C 1.4-1. >*«2*.6*Cl. -021 


GO TO 20 

CC0H1 C 1 1 = C1 *-02)/ CS£*RTCI.-.6*.6>-*-.£s*.6*C .5-02/2 * 
$ Cl, ♦SORT Cl. -•£*•&))) 

23 CC0MlCI)=C2./Ci.4*.6*.c>)*CCRFftC**RGAM>-l.) 

2 CONTINUE 
1 CONTINUE 
RETURN 
END 


OH 1 GH STAL' PAnry m 

0 $ PQ or, J A ^ JST 

)/ ^ Quality 


4 * 44**4 4 4444444444 * 4 * 4444444*4 4*4444 444 4*4444444 4444 * 44 4**4 44 * **** 

SUBROUTINE MATC0HCA,X1,PHU2> 

** #******■**■******■•***■*■**•***■*■***■*'**★*■*•*****•**■**'**'**********'******** 


MATRIX CALCULATION 

DIMENSION PHU1C2C51,X1C2C5) 
DOUBLE PRECISION ACE :»E1)»XK,Z 


POLYNOMIAL CALCULATION 
RHK=-1.*2./21. 

DO 5 IK^l.20 
XK=CBL* (RHK> 

CALL CORR CRHK*PHK ,X1 «PHU2, 2 CO 
CALL DFDXCZ,2.XK,11 
ACIK»1)=C1.D*C> 

DO G J=2,22 

’ A C IK*-21 > = cTl *0*3) ♦DELE CPHK) )*Z*CC*Gu*0) ) *C C.GO-* C )/( 2.D+C ) ) * 

* CC2 .D+-0) ♦DELE CPHKJ 3 * CC 1 .D^ 1 1 ♦DELE CPHK 11* 

% m.o+i)+z*z) 

i RHK=RHK*2./21. 

RETURN 

END 


SUBROUTINE COMPLOT CX1,CC0M1 ,CC0M2, CC0M3 ,CCQM4,X3 ,CP ,N 2 , 
% K AI * ALPHA* IPOS* I POLY * N,R 1 


PLOTTING ROUTINE 

DIMENSION XI C2 - 51 * CC0M1 C 205 1 »CC0H2 C 205) *CCQM3C2C5) , CCQM4C2C5) 
DIMENSION X3C2:5),CPC2r5) 

DOUBLE PRECISION ALPHA, R 

CALL PLOTS 10, 0*4LPLOT> 

CALL FACTOR C. 951 
NI=2LQ 

Xl(Nl+l)=-i. 

X1CN1*2>=.4 

X3CN2+1)=-1. 

X3lN2+2>=.4 
CPCN2+1 )=1. 

CPCN£+2)=-.3 
CCQM1 CN1+1 1 = 2. 

CC0MlCNl-*-2) = CPCN2^2J 
CC0M2 CNl^l ) = CPC N2+1 1 


noon on non non on non 


CC0M2!Nl-*-2>=CPiN2*£> 

CC0H3!Nl*i>=CPiN2+l> 

CCGH3tNlf2>=CPCN2*2> 

CC0M4!Nl«-lJ=CP!N2«-i> 

CC0K4!Ni-*2>=CP!N2+2> 


ORIGINAL PA GE IS 
OF POOR QUALITY 


AXIS 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 


AND ORIGIN 

AXISC2.2.2.2.1HX.-T .Xi! Nl+1> .X1CN1+2 >> 

AXI S! 2*2»2*2»2HCP»2 .B*L *9C • 3.CC0H1 CNI+1) »CCGH 1CN1+2) > 
AXISC2.2* 8.2»2H ,-l .5* *0**999. *1.) 
AXIS!7.2.2*2.0*0.t..9O..999.»l.> 

ORIGIN (2*2 *2*2*0 

LINE CXI *CCOHi*Nl*l »-&»!> 

SYMBOL! 3. .5*5. . l.I *0.»-l> 

SYMB0LC3.4 *5.5* *1* *N=*.0..£ > 

NUMBER C 3*8 *5.5* *I»FLOAT!N> * 3* . 0) 


CALL LINECX1 *CC0M2*N1»1 *-6»2> 
CALL SYMBOL! 3*«5*2E**2.2.3* *“I> 
CALL LINECX1 .CC0M3 .Nl.l *-6 ,3) 
CALL SYMB0L!3*.5»* *1.3.0*.-ii 
GO TO 7 

CALL LINECXl,CC0M4,Nl,l.-6*4> 
CALL SYMBOL! 3* » 4*75* *1*4 * 3 • *— 1> 
CONTINUE 


COMPARISON 

CALL L1NECX3 .CP »N£ .1 . - *£> 


TITLE 

CALL 0 SY^B0L!1...5*.1*»KARHAN-TSIEN CORRECTION FORMULA* . 

S3. ,31) 

21 CONTINUE 

IF !KAI*EQ*£) CALL SYMBOL !. 2 .G.25 »* i4 
*,*NACA C 312 M=D*6 ALPHA=0 **i‘* *25) 

CALL PLOT!-* 5.” *5* ?99> 

RETURN 

END 

,* ****** *********************** ****** ***************** ********* 

SUBROUTINE C0RRCX.UXS.X1.UX.N1> _ . 


CALCULATES UXS BY INTERPOLATION USING UX12G5) 

DIMENSION XI !205>.UX!2 r 5> 

UXS = ! X”X1! 2 > >*! X-XX ! 3> >* CX-X1 14 )>*UX(1)/!X1! 1>-X1 ! 2 > > 

l!X-Xl!i>>*!X-XlC3>>*lX-XlC4>>*UX!2>/tXlC2>-XHl>>/CXH2>-Xl!3>> 

$!X-Xi!l>>*!X-XiC2>>*!X-Xl!4)>*UX!3)/!Xl!3>-Xltl>>/!XH3>-Xl!2>> 

l(X-Xl!" >>*CX-XlC2>>*!X-Xl!2>>*UXt4>/CXlC4>-XHl> >/ ! Xl ! 4 > -XI C2 > ) 
$/! XI ! 4 >— XI ! 3 >> 

ENOIF 


c 


c 


DO 1 I=3»Nl-3 

IF(X*GT.X1CI1 .AND* X.LE.Xi €!♦! > > THEN 

UXS=CX-XlCl))*<X-XlCI*i))*CX-XlCI^2))*UX(I-ll/CCXllI-l)-XHI))* 

*CXiCI-l)-XiCI+i))*CXiCl“i>-Xl(I+2)))^ 

*CX-XiCI-l>>*(X-XlCI*U>*(X-Xi(I+2>>*UXCI)/tlXlCl)-XlCl-l) 1* 
»CX1< I)-XitI^l>)*CXlC I1-X1C I+2J ) )♦ 

StX-XlCl-i>>*<X“Xl<I>>*4X-XlCIf2)>*UXCI*l>/C CX1(I-H>-X1<I-1> >* 
*1X1 ( I+i)-XldJ>*tXlC 1*1 > -X 2 CI+2J a* 

SCX-XICI-i) )* (X-X1CI) >*CX-XiCI«-i))*UXCI-*-2>/t CX1U*2)-X1(I“1>>* 
*CX1 (I+2)-XltI)J*<XlC 

END IF 

1 CONTINUE 


IF(X.GT.X1CN1-Z1> THEN 

UXS = tX-XlCI>>*<X-X2CI«-i> J*<X-X1 C 1*2 I > *UXCI -l)/t C XI C 1-1 >-Xl< I))* 
*<XlCI-l)-XlCI«-i)J*CXltI-l>-XltI«-2J>>* 

*(X-XlCI-l>>*<X-XlCI+l)>*CX-XltI^2>>*UXtI>/CtXlCI)-XlCI-l) )* 

SC XI C I ) -X1C I*-1>>*CX1CI)-X1C I+2) > ) ♦ , , % 

SCX-XiCI-l>>*CX-XlCIJ>*<X-XiCI*2>>*UXCl*l>/CCXlCI*l>-XlCI-l>>* 
SCX1 C I*1J-X1(I))*(X1( I-*-! 1-X1 CI-*-2 > ))♦ 

SCX-X1<I-1>>* CX-X1CIJ >*<X-XlCI+l>)*UXCI«-2>/t CXiCl+2)-XltI-l) )* 
SCX1 CI+2J-X1C I) >*CX1 ( I*2)-XiC I*l»)) 

END IF 


RETURN 

END 


ORIGINAL PAGE IS 

OF POOR QUALITY 
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r* ********************************** * * **★****♦♦*★*♦** ** ** *** *♦★*** **** 

SOLUTION OF THE INVERSE DESIGN PROBLEM 


PROGRAM RE V4 (ANSWER * OUTPUT .COEFPO* CPDAT A*CPTAB4 * JOUKQ 9 PLOT » 
tTAPE2=ANSyER *T APE3=CCEFPG * TAPES =CP DATA .TAPE S=CPT AB 4 *TAPE7=J0UKC > 
DIMENSION XA (5 ) *X2 ( 2 r 5 ) * CP (2 05 ) *UX (2 05 ) *Y1 (20a > » Y2 (2 05 > 


UI M EN§I ON Xfl 13 J * A 2 t 2. ' - * * tr . Jilt 
DIMENSION XDi(205>*VXC2C5).UC2C5> 
DOUBLE PRECISION X *A €5 0 9 51 > .DETER 


INPUT GU ANTITIES 

POLYNOM OF DEGREE fi 
N=24 

PRESSURE DISTRIBUTION 


IfTHEORETICAL ELLIPTIC DISTRIBUTION* ALPHA=i) 

3. POLYNOMIAL NACA0C12 (ALCP* COEFPO N=i9 > * ALP HA=0> 

4, ELLIPTIC POLYNOMIAL (ALCP. COEFPO N=3 0 ) *ALPHA=G> 


ALPHA= : 


8 *CYLINDER*ALPHA=0) 

KA 1=5 

DISTURBANCE FIELD: 1*UX 

2.VX 

KAP=2 


COMPUTES OR READS THE X -DISTURBANCE VELOCITY FIELD 

THEORETICAL ELLIPTIC DISTRIBUTION 

THICKNESS RATIO 

RT=.12 

IF(KAI.EQ.I) THEN 

N 1=2 CO 

XI (1 )=-!• ♦2* /FLOAT12 CD 

UXCI) = ?1* + RT)/ C1.*-RT*RT*X1 CI)*X11 1 >/(l«-Xl(I>*Xl(I)>)-l» 

VX(I>=-(UXCI>*2.>*RT*X1CI>/3GRT(1.-XI(I>*X1(I>> 

Xl(I«-l)=XlCl>-*-2./FLOAT(20‘l) 

END IF 

CYLINDER „ _ 

IF1KAI.EQ.S) THEN 
N 1 = 2 DO 

Xl(i)=-i.«*2./FLOATC2 01> 

DO 30 1=1 »N1 

UX(I)=2.*(1.-X21I)*X1(I)>-1* 

VX(I) = -2**XI (I )*£QRT1 1.-X1 (I)*X1( I>> 

Xl (1+1 )=X1 ( I )^2«/FL0 AT (20 1 > 

END IF 


DATA OF CP FROM PROGRAM: ALCP 


ORIGINAL PAGE IS 
OF POOR QUALITY 


oononno onnnooo 


52 


IF (KAI»EG»3 .OR. KAI.E&.4) THEN 

Nl=20? 

DO 52 1=1 *2.0 

READ (3*255) X1(I),CP(I) 

XS=X2 Cl) 

CALL*DFDX(Z*KAI,XS*RT) 

UX(I) = SQRTC(i.-CPCI))/(I-*-Z*Z) >-i. 

VX(I) = (UX(I)«-1.)*Z 

CONTINUE 

END IF 


ORIGINAL PAGE IS 
OF POOR QUAUT5Q 


FILE OF CP CTAB) AT S DEGREE 
IFCKAI.EQ.5) THEN 
Nl=40 

DO 23 1=1 *N1 

READ (5*255) XI (4i-I),CPC4l-I) 
Xl(41-I)=2.*(Xl(41-I)-.5) 

XS=X1(41-I) 

CALL DFDXCZ*KAI,XS»RT) 

UX (4 1- I) = SORT C (1. -CP (41-1 )>/(!. «-Z*Z))-l. 
23 VX{41-I) = (UX(41-I)-*-l.)*Z 

END IF 

JOUKOUSKY AIRFOIL 
IFCKAI.EQ.7) THEN 

N 2 =2 C C 

READC7* 215) XK1),CP(1> 

DO 57 1=1 *25 G 

READ (7*2 25) Xl(I).CPd) 

XS=X1(I> 

CALL DFDX(Z*KAI,XS,RT) 

UX (I ) = SQR T( ( 1 • -CF ( I) )/ Cl •♦Z*Z ) ) “1 • 

VX(I) = (UX(I)-*1.)*Z 
57 CONTINUE 

ENDIF 


** ******** ■a****-********************************* **** * *************** 


COMPUTES THE MATRIX 

X=-(l.D«-3) + C2.D*Q)/DBLE(FL0AT(N«-l)) 
DO 26 1=1* N 
XS=SNGL(X) 

CALL CDRR (XS *VXS* XI* VX* Nl) 

CALL CORR (XS »UXS»X1*UX*N1) 

CALL MATR1(A*I « X* UXS * VXS*N*K AP) 
X=X*(2.DfG)/DeL£(FLCAT(N*l)) 

26 CONTINUE 




** ****»•»**♦■***■*****■*•*■***-*+** *********** ♦*♦■*♦**•**■****■* *********** * * * 


SOLVES FOR THE MATRIX A(N.NPLUSM) 

CALL GAUSS C A, N. DETER) 

WRITE (2*223) DETER tN »N*1 
WRITEC2* 22e) 

00 5 1=1 *N 

XA ( I)=SNGL( A (I *N*1) ) 


ooonnn o nonnn n on oooonoo 


. .. 


UR ITE C2 * 2 31) XA (I } 
5 CONTINUE. 


. * O************* A***************** 


A********************* ******** 


COMPUTES THE AIRFOIL PROFILE USING THE POLYNOMIAL COEFFICIENTS 
PI=3. 141592654 

IFCKAI.EQ.2 .OR • KAI.EQ.3 .OR. KAI.EQ.5 .OR. KAI.EQ.6 
% .OR. KAI.EC.7) THEN 
N0=2S1 
Y1C1>=C. 

X01 Cl )=1. 

00 28 I=2»2C1 

XDltI>=XDlCl-l>-2./FL0ATC25i> 

XS=XDltI> 

IFCKAP.EQ.2 > CALL C0RRCXS,VXS*X1 »VX»N1> 

IFtKAP.EQ.l) CALL C0RRCXS*UXS*X1 *UX»NI> 

CALL CPDRlCXS»XA»UXS»VXS*N.KAP> 

FX=VXS/tl*-K»XS> 

Y1CII=Y1CI-1>^FX*CXC1CIJ-XD1 CI-D) 

28 CONTINUE 
END IF 

ELLIPTIC AIRFOIL _ . 

IF C KAI «EQ. 1 .OR* KAI *_G. 4 .OR. KAI.EQ.S) THwN 

nd=i*:c 

XD1 C 1>=-1./FL0AT(2 '1 > 

Y1C1>=.12 

IFCKAI.EQ.fi> Y2C1>=1. 

DO £•'; 1=1*99 

XD1C I-*-l>=XDl CI> “2 ./FLOAT (221 ) 

X 3 ~XD 2 ( ^ J 

IFCKAP.lQ.II CALL CORRC XS,UXS *XI *UX* Nl> 

IF (KAP.EQ.2) CALL CGRRtXS»\IXS »X1 *VX»N1> 

CALL CPdR1(XS*XA*UXS*VX3*N*KAP> 

FX=VXS/Cl.+UXS> 

YlCI*l>=YitI>-»-FX*CXCl tl+D-XDIC I>> 

£ CONTINUE 

END IF 


4 

COMPUTATION OF THE THECRITICAL PROFILE FOR COMPARISON 


CALL THE0CND*KAI»X01*Y2#RT> 


2 C2 
2 S3 
2 Cl 
2 08 
2 25 
216 


FORMATC/ *• SHALL PIVOT - MATRIX MAY BE 
FORMAT C/** DETER = **C16.12/»*N = 

FORMAT ( 1H* 1 3F1C .5) 

FORMAT C/*POLYNOMIAL COEFFICIENTS •/> 

F0RMATC2F1..6J 

FORMATCI3) 


SINGULAR*) 
*.I2/»*NPLUSM = 


• * 12/) 




PLOT OF THE PROFILE OF THE WING 


ORIGINAL ? AGS IS 
OF, POOR QUA.UTY 


| 


non n n n o nonnnon 


CALL PLOTS !C ,0,4LPLO7) 

CALL FACTORS. 95) 

CALL SCALE ! XC1 ,5. ,ND *1) 
XDi!ND+l)=-l. 

XD1 !NQ+2)=«4 

IFCKAI.EQ.l .OR. KAI.EO.4 .OR. KAI 
CALL SCALE!Y1,6.,ND,1) 

Y1!ND*I)=0.0 
Y1!ND*2>=. ?3 

IFIKAI.EQ.8) Y1 (NO* 2 )=.2 
Y2!ND*l)=Yl!ND*i) 

Y2!ND*2)=YHND*2) 

CALL AXIS ( 2 .2,2.2, 1HX»-1 ,5»*G»,XD2 
CALL AXIS!2.2,2.2,2HY ,2,6.,S0.,Yi 
CALL AXIS!2.2,8.2,:,C,5.,C.,999.,I 
CALL AXIS!7.2*2.2,0,O»£.,9C.,999., 
CALL 0RIGIN!2.2,2.2,D) 

CALL LINE! XD1, Y1,N0, 1,-6 ,1 ) 

CALL LIN£!XD1,Y2,ND, 1,0,2) 

CALL SYMBOL ( .6,5.5, .1,1, 0. ,-1 ) 

CALL SYMBOL !1«,5»5,«1, 'POLYNOMIAL 
CALL NUMBER !2.5,5.5,.i,FLOAT!N),3. 

TF (Kir.rD.' ! i rill c YHFtni i _ r? = 


ORJGliJ* T' r r 

OE POOR Q C/ALi 


•EG.8) XD1!ND*2)=.2 


CALL AXISC2. 
CALL AXISC2. 
CALL AXIS(7 • 
CALL ORIGIN! 
CALL LINECXO 
CALL LINEiXD 
CALL SYMBOL! 
CALL SYMBOL! 
CALL NUMBER! 
IFtKAI.EQ.2) 

53.. 26) 

IF !KAI.EG.4) 

53., 26) 

IF !KAI.EQ.E) 

50.. 24 ) 

IF !KAI.Eq.1 ) 
IF ! KAI .EQ.8) 


!ND*1),XDUN0*2)) 

!ND*1),Y1!ND*2)) 

. ) 

1.) 


CALL SYMBOL!, 


• 6.25. 


N=»,0.,12) 

, C) 

14, *DATA FROM THE PROGRAM 


IF !KAI.EQ.7) CALL SYMB 
CALL PLOT!-. 5,-.5, 999) 
CONTINUE 
ST CP 

END 


CALL SYMBOL!. 2, 6.25*. 14, 'DATA FROM THE PROGRAM ALCP*, 

CALL SYpeCL!. 2 , 6 . 25 ,.i 4 ,*DATA FROM TAB !3 DEGREE)*, 

CALL SYMBOL! *2 ,6.25 , .14, * ELLIPTIC AIRFOIL *,0., 16) 

CALL SYMBOL t .2 ,6.25, .14, • CYLINDER * , C. ,8 ) 

CALL SYMBOL! .2 ,6.25 ,.14, * JOUKOUSKY AIRFOIL** 0. ,17) 


14, 'ELLIPTIC AIRFOIL*, 0.,: 
14, • CYLINDER * , C. ,8) 

14, 'JOUKOUSKY AIRFOIL*, 3. 


CALCULATES THE AIRFOIL DERIVATIVE 

SUBROUTINE DFDX!Z,KAI,X,RT) 

ELLIPTIC AIRFOIL 

IFIKAI.EQ.1 .OR. KAI .EC. 4) Z=-RT*X/ !SQRT !1.-X*X) ) 

CYLINDER J 

IF(KAI.EQ.B) Z=-X/SQRT!1.-X*X) 

NACA ^012 

IF (KAI.EQ.3 »CR. KAI.EG.5 .OR. KAI. EG. 6) THEN 
Z=1 7 • +RT* ! .174225/! SORT ! !X*1 • )/2.) )-• 063 
S -. 3516*!! X*l.)/2.)*. 42645* !!!X+l.)/2.)**2.) 

S -.273*!!! X* l.)/2»)**3.)) 

END IF 

JOUKOUSKY AIRFOIL 

IFIKAI.EQ.7) Z=-4.*RT*!SGRT!i.-X*X)*!i.-X)*X/SGRT!l.-X*X) )/3./ 
SSGRT !3. ) 

RETURN 

END 


no on nnnnnnn ono 


CALCULATES UxS BT INTERPOLATION USING UXC2051 


SUBROUTINE CORR (X.UXS.Xi «UX»N1> 

DIMENSION XiC2G51tUX€2C51 
IFCX.LE.X1C311 THEN 
UXS=CX-Xlt21 1*CX“X1C31>* IX-X1C4 1 1*UXC11/<XI 

S/CX1(1)-X1C3)>/CX1C 1>-X1CA))> 

SCX-XKil l*<X-XlC31T*CX-XlC4>l*UXC2>/(XiC21- 
S/CX1C21-X1C411* 

SCX -XI (ll)*CX-Xlt21 1*CX“X1 C 4 1 1*UX (31 /(XI (3) - 
$/(Xl(31rXl(4 ))♦ 

SCX-X1 III 1*( X“X1 (21 1* (X“X1 (3 1 1*UX (41/1X1(41“ 
S/CXH41-X1C31) 

END IF 

DO 1 I=3»Nl-3 

IF(X.GT.XI(I1 .AND. X.LE .XI 11*1 1 1 THEN 
UXS=CX“XlCl>)*CX-XlCI^l>)*tX“XiCI'*-2)l*UXCI- 
S(Xl(I-ll-Xl(I«-lll*(XltI“ll“Xl(I>»-211!* 

S(X-X1 (1-111* (X-Xl(I*lll*(X-Xl(I«-211*UX(Il/( 
$(X1(I1-X1(I+111*(X2( I) “XI ( I*21 1 1 ♦ 
SCX-Xl(I-l) 1* (X-XKIl 1*<X“X1(I«>2} 1*UX( 1*1 1/( 
$(Xi (I*1)-X1( I) 1*(X2( I«-2 1-X1 «!♦£)>>♦ 
SCX-XlCI-i)>*CX-XltU 1*(X-X1(I«*' 11*UX( 1*2 1/( 
S(X1 C I+21-XK 111* exit I*2)“XlCI«-2>>> 

END IF 
l CONTINUE 

IF(X.GT.X1(N1-211 THEN 
I=Nl-2 

UXS=(X-X1CI1 l*CX“Xi€ I^111*<X-X1(I-*'21}*UXCI“ 
sexicl-ll-xli I«-111*(X2CI-11-X1(1-*2111* 
S(X-X1(I-111* (X-XI(H>lll*(X-Xl(I*2>l*UX(I>/( 
S(XltIl-XlCI«’lll*CXlCIl“XICIf2>l)* 
$(X-X1(I-111*(X“XKI11*CX-X* (I«-211*UX(I*ll/( 
scxici^i}“xiciii*cx:ci^ii-xiti^2iii^ 
S(X-Xl(I-lll* CX-X1CI1 }*(X-Xl(I*ill*UX(I*21/( 
*lXltI^21-XHIll*CXItI+2l-XlCI>llll 
ENDIF 
RETURN 
END 


1 1) -Xi C 211 

Xl(lll/(X1(21-Xl(311 
X1C11 )/ ( XI C 31“XI C2 1 1 
X1C11 1/CXH41-X1C211 

11/CCXlCI-ll-XlCIl)* 
(XKIl-XKI-ll 1* 
(X1(I«-11-X1(I-11 1* 
(XKI+21-XKI-11 1* 


11/((X1(I-11-X1(I11* 
(XKIl-XKI-lll* 
(X1(I*11-X1(I-11 1* 
CX1CI+21-X1CI-1) 1* 


** ****** ************ **************** ****************** * ************* 


CALCULATES THE THEGRITICAL PROFILE OF THE AIRFOIL 

SUBROUTINE THEQ (Nl»KAItXl*Y2*RTl 
DIMENSION X1C2C51*Y2(2:51 

ELLIPTIC AIRFOIL 
IF ( KAI .EQ. 2 .OR. KAI.EC.41 THEN 
DO I I=i»Nl 

1 Y2(Il=RT*SGRT(l.-Xi(Il*Xi(Ill 

ENDIF 

CYLINDER 

IF (KAI.EQ.9 1 THEN 
DO 10 1=1 *N1 

1? Y2(Il=SQRT(l.“Xi(Il*Xl(Ill 

ENDIF 


on on n o n o o nnnnnnnnnn on 


C 


NACA "012 

IF (KAI .EG.3 .OR. KAI.ES.S -OR. KAI.EQ.b) THEN 
DO 2 I=1»N1 

Y2CI>=i:.*RTM.2969*S£lRTCXlCI>*.5-»-.5>-.i26*CXlM>*.5«-.5J 
C -•3516*CCXItn*.5*.EI**2)«-.2843*CtX2CI)*.5*.5>**3) 

C -.i:’5?*((Xi(I)*.5+.5)**4)) 

2 CONTINUE 
ENDIF 


4 


JOUKUSKY AIRFOIL 
IF CKAI.EQ.7) THEN 
DO 4 1=1 »N1 


BKIGmax; page 
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Y2CI)=4.*RT*Cl.-XlCI>>*CSQRT(l.-XICI>*XlCI)>)/3./SQRT<3.) 


CONTINUE 


ENDIF 


RETURN 

END 


** «*•#*** **********♦•.**********•********. ******■*♦•.■**•** **».** ********* 


GAUSS -JORDAN REDUCTION 

THIS SUBPROGRAM FINDS THE SOLUTION VECTOR CORRESPONDING TO A 
SET OF N SIMULTANEOUS LINEAR EQUATIONS USING THE GAUSS-JORCAN 
REDUCTION ALGORITHM UITH THE DIAGONAL PIVOT STRATEGY* 

SUBROUTINE GAUSS (A .N .DETER ) 

DOUBLE PRECISION ACS '.51> .DETER .EPS 
EPS=l.D-i: 

NPLUSM=N-*i 

BEGIN ELIMINATION PROCEDURE 

DETER=Ci.D+G> 

DO 9 K=1.N 

UPDATE THE DETERMINANT VALUE 

DETER=DETER* AtK.K) 

••••• CHECK FOR PIVOT ELEMENT TOO SMALL .*••• 

IF(DABSCACK.K)) *GT *EPS> GO TO 5 
GO TO 7 

••••• NORMALIZE THE PIVOT ROy ..... 

5 KPi=K*l 

DC fs J=KP1 »NPLUSM 

6 A(K.J>=A(K.J>/ACK.K) 

ACK.K>=1.D«-C 

••••• ELIMINATE KCTH) COLUMN ELEMENTS EXECPT FOR PIVOT ••••• 
DO 11 1=1. N , * 

IF < I.EG.K .OR. ACI»K).EG.(C.D«-0) > GO TO 11 
DO 8 J=KP1 .NPLUSM 

8 AC I. J)=A(I. J 1-A C I » K) *ACK« J ) 

ACI.K>=0.D*C 

31 CONTINUE 

9 CONTINUE 
GO TO 1 

7 CONTINUE 
STOP 

1 CONTINUE 
RETURN 
END 

** ****************** ********************************************** 
COMPUTES THE MATRIX ELEMENTS 


K ba 


non o o no 


SUBROUTINE MATR1CA*I *X»UXS *VXS» N»H«P) 
DOUBLE PRECISION X.A <5€ *51 > *PI 
PI =3 .141 5926 54D+0 


IFCKAP.EQ.l) THEN 

A ( I * 1 )=CDL0& C C 1 1*D* - J ♦Xl/i C i .D+C )-X> > )/PI 
DO 2. J=2 t N 

ACI,J)=X*ACI*J-'> -CC1.D+C3>-DBLE CFLOATC C-l)** CJ-1>) >>/ 
* COBLE CFLO AT CJ-1 > ) *PI> 

1 continue 

ACI,N+1)=DELECUXS> 

ENDIF 


IFCKAP.EQ.2) THEN 
AC 1*1 >=1.0+0 
DO 2 J=2*N 
ACI»J>=X**CJ-1> 
CONTINUE 

AC I * N+l) =DBLE C VXS) 

ENDIF 

RETURN 

END 


ORIGINAL PAGE IS 
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PRESSURE DISTRIBUTION 


C 


C 


c . 


SUBROUTINE CPDRlCXS,XA,UXS,VXS t N,KAP> 
DIMENSION XAC52> 

IF CKAP.EQ. 1) THEN 

VXS=XAC1J 

DO 1 1=2 *N 

VX S=VXS+XACI>*CXS**CI-i) > 

1 CONTINUE 
ENDIF 


IF CKAP.EQ. E> THEN 
PI =3 .1415920540+0 
RINT=AL0GCC1.+XS)/C2 .-XS))/PI 
UXS=XAC1>*RINT 
DO 2 J— 2*N 

F, INT =XS*RINT-C 1. -FLQATt C- 1 >+*CJ~l>)> / C FLO AT l J“l> ) /P I 
uxs=uxs+xacj>*rint 
CONTINUE 

ENDIF * 


RETURN 

END 


